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ABSTRACT 


A  Fourier  transform  technique  for  inverting  gravity 
data  has  been  developed  to  yield  a  two  dimensional 
multilayer  density  profile.  The  method  extends  an  algorithm 
originally  developed  by  Parker  (1972)  and  Oldenburg  (1974). 
This  multilayer  inversion  scheme  requires  an  initial 
starting  model  from  a  seismically  determined  model  and 
allows  for  lateral  variation  of  density  and  thickness  within 
each  stratum.  The  inversion  procedure  is  used  to  study  the 
crustal  structure  of  the  southern  plains  of  Western  Canada 
and  the  resulting  inversion  models  are  in  reasonable 
agreement  with  additional  independant  seismic  models.  It  is 
shown  that  the  algorithm  may  be  rearranged  to  find  the 
density  distribution  of  a  near  surface  slab  in  order  to  fit 
the  residual  gravity  anomaly.  This  appears  to  be  a  valuable 
tool  for  interpreting  shallow  features  responsible  for  the 
short  wavelength  components  in  the  Fourier  transform  of  a 
gravity  profile.  The  density  distribution  at  the  top  of  the 
Precambrian  crustal  basement  across  the  Churchill^Super ior 
provinces  has  been  computed  and  it  appears  to  be  useful  in 
interpreting  the  position  of  the  boundary. 
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CHAPTER  1  INTRODUCTION 


1.1.  Previous  work: 

Many  methods  have  been  developed  to  determine  the 
shape,  density  distribution  and  depth  of  a  subsurface  body 
from  a  measured  gravity  anomaly.  In  the  thirties  and  forties 
with  the  introduction  of  use  of  graticules,  dot  charts  and 
other  similar  graphical  aids  to  compute  gravitational 
attraction,  attempts  were  made  to  perturb  an  initial 
starting  structure  until  the  gravity  calculated  matched  the 
observed  field.  (Skeels,  1947)  In  theory  these  methods  can 
be  made  as  precise  as  one  pleases,  merely  by  increasing  the 
scale  to  which  the  graticule  is  constructed,  however,  in 
actual  practice  this  may  be  difficult,  if  not  impossible. 

Talwani  et  al  (1959)  devised  the  line-integral  method 
to  compute  the  gravity  of  a  perturbing  body.  Talwani's 
inversion  scheme  involves  approximating  the  cross-sectional 
shape  of  the  perturbing  body  by  polygons.  The  residuals 
between  the  calculated  and  observed  fields  are  used  to 
adjust  the  model  parameters:  density  contrast,  shape,  and 
depth.  By  choosing  a  sufficient  large  number  of  sides  for 
each  polygon,  the  cross-sectional  shape  of  the  body  can  be 
approximated  to  any  desired  accuracy.  The  main  drawback  of 
this  method  is  that  the  interpreter  has  to  decide  what 
changes  to  make  for  the  shape  of  the  model.  As  the  shape  of 
most  perturbing  masses  is  not  known,  such  changes  will  bias 
the  final  interpretation. 
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Veiling  Meinesz  and  others  (1934)  proposed  the  use  of 
rectangular  blocks  of  variable  size  for  the  computation  of 
gravity  of  any  irregularly  shaped  bodies.  Bott  (I960) , 
Corbato  (1965),  Tanner  (1  967),  Negi  and  Garde  (1969) 
approximated  the  initial  model  by  a  set  of  rectangular 
prisms  of  constant  density,  and  calculated  its  gravitational 
field.  The  differences  between  the  calculated  and  observed 
fields  are  used  to  readjust  the  heights  of  these  rectangles 
(Corbato,  1965) .  This  method  can  be  made  as  precise  as  one 
pleases  by  using  a  sufficiently  large  number  of  small 
blocks.  However,  as  the  number  of  blocks  is  increased  the 
computations  become  increasingly  tedious,  and  the  iteration 
scheme  become  unstable  and  converge  slowly. 

Parker  (1972)  derived  an  algorithm  to  compute 
gravitational  attraction  of  a  three-dimensional  source  using 
Fourier  transform  method.  The  use  of  the  Fast  Fourier 
Transform  or  FFT  (special  issue  of  IEEE,  1967)  resulted  in  a 
saving  in  computation  time  compared  with  other  methods  in 
the  analysis  of  three-dimensional  gravity  sources.  Parker's 
three-dimensional  scheme  can  be  simplified  to  obtain  two- 
dimensional  gravity  profiles. 

Oldenburg  (1974)  rearranged  Parker's  two-dimensional 
formula  to  obtain  an  iterative  procedure  to  calculate  the 
shape  of  the  surface  between  two  constant  density  media. 

The  purpose  of  this  thesis  is  to  expand  the  Parker- 
Oldenburg  inversion  algorithm  to  many  layers,  allowing 
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densities  to  vary  along  the  length  of  the  gravity  profile. 

The  inversion  of  two  gravity  profiles  in  the  southern 
parts  of  Manitoba  and  Saskatchewan  is  offered  as  a  practical 
application  of  this  multi-layer  inversion  scheme  (Fig  3.1). 

This  modified  version  of  the  Parker-Oldenburg  inversion 
algorithm  is  later  rearranged  to  find  the  density  variation 
of  a  near  surface  slab  in  order  to  fit  the  difference 
between  the  original  gravity  and  the  gravity  computed  from 
the  inverted  topography. 

1.2.  Parker-Oldenburg  Algorithm 

The  derivation  of  Parker*s  (1972)  formula  was  fully 
discussed  in  his  original  paper  and  subsequently  by  Hinson 
(1976).  A  full  discussion  of  the  derivation  is  presented  in 
Appendix  A. 

Parker's  algorithm  for  a  three-dimensional  source  with 
constant  density  is: 

^  00  _ 

FCAg]  =  -2WSd  exp  (- 1  K |  ZQ)  X  (ill  n-»/n!)  F[  hn  (r)  J  (1) 

n=  i 

where  F[  A  9]  is  the  Fourier  transform  of  the  gravitational 
anomaly,  G  is  the  universal  gravitational  constant,  TT  is  the 
wavenumber  vector  on  the  x-y  plane,  T  is  the  position  vector 
on  the  x-y  plane  from  the  origin  to  the  source  element,  d  is 
the  density  ,  Zo  is  the  distance  from  the  origin  to  the 
observation  datum  plane,  and  h(r)  is  the  upper  boundary  of 


. 
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the  perturbing  mass  measured  from  the  x-y  plane.  (See  Figure 

1.1) 

For  a  two-dimensional  source  the  gravity  A  g  is  now  a 
function  of  x  not  IT,  and  TT  is  now  a  scalar  not  a  vector.  To 
obtain  Oldenburg's  inversion  scheme  the  first  term  of  the 
infinite  sum  is  transposed  and  rearranged. 

© 

F[  h  (x)  ]  =  -exp  (KZo)  F[Ag  (x)  ]/27TGd-  T  jK^/n!)  F[hn  (x)  ]  (2) 

n~2 

Assuming  Zo  and  d  are  known,  andAg(x)  is  given,  this 
expression  can  be  used  iteratively  to  obtain  the  topography 
h  (x)  of  the  perturbing  body.  An  initial  guess  of  h (x)  (even 
for  h(x)=0  is  satisfactory)  is  input  to  the  right  hand  side 
of  the  equation  (2) •  The  infinite  sum  is  computed  until 
convergence  criterion  is  met.  An  updated  value  of  h(x)  is 
obtained  by  taking  the  inverse  Fourier  transform  of  the 
right  hand  side  of  equation  (2) .  This  updated  version  of 

h  (X)  is  then  put  back  in  the  right  hand  side  of  equation  (2) 
and  the  iterative  procedure  is  continued  until  some 

convergence  requirement  is  met  or  a  maximum  number  of 

iterations  has  been  performed. 

The  conditions  necessary  for  convergence  and  the 
convergence  criteria  of  the  infinite  series  and  the 
iterative  procedure  were  both  discussed  in  detail  by  Parker 
(1972)  and  Oldenburg  (1974).  They  are  not  going  to  be 
repeated  here.  However,  a  summary  is  presented  in 

appendix  B. 
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Figure  1.1:  Subsurface  configuration. 
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The  convergence  criteria  of  the  infinite  sum  of 
Oldenburg  (1974)  algorithm  was  Sn/S2<5E-3  while  that  of 
Parker  (1972)  algorithm  was  Rn/R1<10E-6,  where: 


Sn  =  max  over  all  K  |  ( |  K |  n~l/n ! )  F[  hn  (x)  ]|  (3) 

and 

Rn  =  max  over  all  K  1  exp  (- |  K  j  Zo)  ( j  K  |  n”  4/n! )  F£  hn  (x)  J  |  (4) 


As  for  the  iterative  method  the  convergence  criterion 
required  that  the  mean- root-squared  difference  between  two 
consecutive  approximations  of  h(x)  be  less  than  0.5  metre  or 
ten  iterations  were  run.  The  mean-root-squared  error  can  be 
expressed  mathematically  as 


Error 


N 


Z  [  hj  (x)-hi+l  (x)  Y 


i  =  1 


(5) 


Parker  has  shown  that  for  Zo>0  the  rate  of  convergence 
of  the  algorithm  is  most  rapid  if  H/Zo  is  a  minimum;  where  H 
is  the  maximum  value  of  |h(x)|.  However,  if  H/Zo  approaches 
1  or  greater  than  one  the  algorithm  may  still  converge,  but 
at  a  slow  rate.  To  make  H/Zo  minimum,  the  origin  plane 
should  be  placed  at  the  median  of  h(x).  Since  the  initial 
placement  of  the  origin  plane  is  entirely  arbitrary,  a 
displacement  of  the  origin  plane  does  not  affect  the 
validity  of  equation  (1),  but  it  does  alter  the  numerical 
values  of  Zo  and  h(x)  and  thereby  H. 


1.3.  Filtering 


It  was  indicated  that  straight  forward  application  of 
the  iterative  procedure  usually  resulted  in  divergence, 
especially  when  the  relief  of  h  (x)  was  high.  To  solve  this 
problem,  the  standard  technique  of  filtering  downward 
continued  data  was  applied  to  the  iterative  procedure.  Since 
most  short  wavelength  anomalies  are  caused  by  shallow 
features,  it's  therefore  justified  to  remove  the  high 
frequency  components.  To  remove  these  high  wavenumber 
oscillations  the  right  hand  side  of  equation  (2)  is 
multiplied  by  a  low  pass  filter  B (K)  which  cuts  off  all 
spatial  frequencies  above  k(WH)  and  passes  all  spatial 
frequencies  up  to  k(SH)  .  The  filter  used  by  Oldenburg 
(1974)  was: 

1  ;  |  K/27rj  <  k  (WH) 

B  (K)  =  1/2  {1  +Cos[  (K-27Tk  (WH)  ) /2  (k  (SH)  -k  (WH)  )  ]} 

;k  (WH)  <1  K/2^!<k  (SH) 

0  ;  |  K/21TJ  >  k  (SH)  (6) 

Wavenumber  K  is  related  to  spatial  frequency  by:  K=27rk. 
This  filter  was  chosen  over  several  other  possible  filters 
because  it  reduces  Gibb's  phenomenon. 


It  was 

also  indicated  that 

judicious 

choice 

of 

k  (WH)  and  k(SH)  could 

be  very 

critical. 

For 

an  y 

gravitational 

profile 

the  filter 

parameters 

k  (WH) 

and 

k  (SH)  can  be  adjusted  to  ensure  convergence  of  the  iterative 
process.  A  decrease  in  the  value  of  k(WH)  and  k  (SH)  will 
result  in  faster  convergence.  However,  if  k(WH)  and 
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k  (SH)  are  too  low,  the  filtering  effect  will  be  severe  and 
the  inverted  topography  will  lose  most  of  its  detail.  In  the 
extreme  case  the  gravity  anomaly  computed  from  the  inverted 
model  might  not  necessarily  agree  with  the  observed  gravity. 
On  the  other  hand  ,  if  k(WH)  and  k(SH)  are  too  high, 
divergence  of  the  iterative  process  might  occur. 

Judicious  choice  of  Zo  and  d  is  vital  to  convergence  as 
well.  If  Zo  is  fixed,  the  number  of  iterations  will  increase 
when  d  decreaes.  Likewise,  if  d  is  fixed  increasing  Zo  will 
require  more  iterations  to  achieve  convergence.  Actually 
there  exists  a  minimum  value,  d(min.),  and  a  maximum  value, 
Z(max.),  beyond  which  the  iterative  procedure  no  longer 
converges. 


•  |  H  v  %  M  :  J 


CHAPTER  2  MULTILAYER  GRAVITY  INVERSION  METHOD 


2.1.  Extension  of  Oldenburg  Inversion  Algorithm 

Oldenburg’s  inversion  scheme  is  based  on  the  assumption 
that  the  observed  gravitational  anomaly  is  caused  by  a 
single  surface  between  two  constant  density  media.  Such  an 
inversion  solely  attributes  the  variation  of  the  gravity 
attraction  to  the  change  of  topography  of  one  interface 
only.  Examples  that  will  satisfy  these  conditions  are  rare. 

As  most  geological  examples  consist  of  many  layers,  any 
variation  of  shape,  bedthickness ,  density,  etc.,  of  beds 
overlying  the  inversion  boundary  can  explain  the  same 
gravitational  anomaly  as  that  explained  by  the  shape  of  a 
single  surface  alone.  Thus  the  inverted  topography  of  a 
given  profile  obtained  by  attributing  all  the  variation  in 
mass  completely  to  a  single  interface  between  two  constant 
density  media  will  be  very  different  from  the  topography 
that  takes  into  account  of  the  effect  of  the  overlying 
layers. 

The  problem  we  are  going  to  study  is  as  follows:  Given 
a  two-dimensional  gravity  anomaly  resulting  from  a  multi¬ 
layered  medium  and  assuming  that  the  thickness  and  density 
of  each  layer,  obtained  from  seismic  information  at  a 
limited  number  of  points  along  the  profiles,  varies  in  a 
known  way  or  else  linearly  along  the  profile  and  that  each 
interface  takes  on  the  same  shape  as  the  bottom-most 
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interface  (called  the  inversion  interface  or  boundary 
throughout  this  thesis) ,  what  is  the  shape  of  the  inversion 
boundary? 

To  extend  Oldenburg's  inversion  procedure  to  allow  for 
variation  of  bedthickness  and  density  of  beds  overlying  an 
inversion  interface,  we  revert  to  Parker's  (1972)  original 
theory.  For  a  three-dimensional  source  with  variable  density 
Parker's  (1972)  equation  is: 

00  _ 

F(Ag  J  =  -27TGexp  (-  |  K|  Zo)  Y  UK|rt”l/n!)F[d(r)hn  (r)  |  (7) 

n=i 

If  the  perturbing  body  is  two-dimensional,  equation  (7) 
becomes : 

00 

F£Ag]=  =-27TG  exp  (-  I  K|  Zo  )  £  ( j  K |  H-i/n! )  F[  d  (x)  hn  (x)  J  (8) 

ns  i 

To  compute  the  total  gravitational  effect  of  a  medium  of 
many  layers  equation  (8)  takes  the  following  form: 


(  Zo  is  replaced  by  Z  for  the  rest  of  this  thesis) 


F[  g  (total)  J 


j=M  00 

-27ns£exp(-|K|Zj)  £  (|K|  n-i/n!)  F[dj  (x)  hj  (x)  j 
j=l  n=i 

(9) 


M  is  the  number  of  layers  in  the  medium.  The  gravity  of  the 
bottom-most  layer  (inversion  interface)  is  designated 
as  g(inv)  ,  and  that  of  the  overlying  layers 
as  g(sum)  .  g  (total)  is  the  summation  of  g(sum)  and  g(inv)  . 


. 


To  evaluate  g  (inv)  M  is  set  equal  to  1  in  equation  (9),  that 


is: 


F[Ag(inv)]  =  -2WG  exp(-|K|Z|)  ^(|K|n-»/n!)F[d,  (x)h"  (x)  ] 

0=1 

(10) 

Transposing  the  n=1  term  from  the  infinite  sum  and 
rearranging, 

F[d,  (x)h,  (x)  J  =  *  (F[£g  (inv)  ]/2TTG)  ex  p  ( \  K  |  Z ,  )- 

GO 

£  (I  K(n“  Vn!)  F[  d,  (x)  h,  (X)  J 

n=2  (11) 

Since£g(inv)  =  g  (total)  -  g  (sum)  , 

F[  d  ,  (x)  h,  (x)  J  =  -F  {[  g  (total)  -  g  (sum)  J/27TG}  exp ( j  K  j  Z ,  )  - 

00 

T  ( I  K|n-1/n! )  F[  d|  (x)  h*J  (x)  ] 

0=2  (11  A) 

where 

j=M  CO 

F[  g  (sum)  J  =  -27TG^  exp  (- |KJ  Z: )  T  (1  K  j^-i/n!)  F[d,  (x)  hj  (x)  J. 

j=2 


Equation  (11)  and  (11A)  is  the  multilayer  inversion 
scheme. 


Essentially  equation  (11)  is  similar  to  equation  (2) 
except  the  densities  at  various  locations  are  allowed  to 


vary. 


Knowing  the  densities  at  various  locations,  the  most 
recent  determination  of  h(x)  is  multiplied  by  d (x)  and  the 
infinite  sum  is  computed  until  convergence  is  met.  Division 
of  the  inverse  Fourier  transform  of  the  right  hand  side  of 
equation  (11)  by  d(x)  will  generate  an  updated  version  of 
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h  (x) .  This  revised  value  for  the  topography  is  substituted 
back  into  equation  (11)  and  the  whole  iterative  procedure  is 
continued  in  the  same  manner  as  previously  described. 

It  is  most  unlikely  that  densities  at  all  locations 
along  the  gravity  profile  are  known.  However,  if  densities 
at  two  or  more  positions  are  known,  either  from  seismology 
or  from  borehole  data,  and  if  densities  can  be  assumed  to 


vary  linearly. 

or  as  some 

known 

functions. 

a  series 

of 

digitised 

densitites 

can 

be 

produced  to 

satisfy 

the 

requirement 

of 

equation 

(11) 

.  To 

remove  the 

gravitational 

effect  of  beds  overlying  the  inversion  surface  (  g(sum))  to 
obtain  a  reduced  gravitational  anomaly  (  g  (inv) )  for 
inversion  we  have  to  strip  off  the  gravity  of  the  overlying 
layers.  To  do  this  the  variations  of  bedthickness,  density, 
depth  and  shape  of  these  beds  have  to  be  known.  Information 
regarding  the  shape  of  these  beds  is  difficult  to  obtain 
unless  detailed  seismic  surveying  has  been  done.  However, 
the  topography  of  these  beds  can  be  assumed  to  change  in  the 
same  manner  as  the  boundary  we  want  to  invert.  This 
assumption  is  justified  on  the  basis  that  the  earth  is 
radially  homogeneous  to  a  first  approximation.  In  the 
absence  of  other  pertinent  information  it  is  the  best 
assumption  we  can  make.  Data  concerning  the  variation  of 
bedthickness  of  each  bed  at  all  locations  is  not  often 
available  either.  However,  often  seismic  survey  or  borehole 
data  will  furnish  enough  information  about  the  average 
bedthickness  and  density  of  each  bed  at  two  or  more 


. 
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locations  to  allow  for  linear  interpolation.  Knowing  the 
average  bedthickness  and  density  of  different  layers,  and 
assuming  the  topography  of  overlying  layers  to  vary  in  the 
same  manner  as  the  inversion  interface  we  can  now  proceed  to 
strip  off  the  gravitational  effect  of  all  overlying  beds, 
and  do  the  inversion. 

2.2.  The  Multilayer  Inversion  Scheme 

The  inversion  procedure  is  illustrated  schematically  by 
the  flowchart  in  Figure  2.1.  The  flowchart  itself  is  self 
explanatory,  however,  it  is  necessary  to  elaborate  on  a  few 
points. 

Only  a  fraction  of  the  observed  gravity  is  used  for  the 
first  inversion  because  the  inversion  horizon  only  accounts 
for  a  fraction  of  the  observed  gravity.  The  choice  of  this 
fraction  is  critical  to  the  number  of  iterations  needed  to 
achieve  convergence.  If  it  is  too  large  or  too  low,  more 
iterations  are  required  to  reach  convergence.  A  good 
estimation  of  this  fraction  can  be  obtained  by  dividing  the 
density  of  the  perturbing  mass  by  the  summation  of  the 
densities  of  all  layers. 

To  calculate  the  gravitational  anomaly  of  beds 
overlying  the  inversion  boundary  it  is  necessary  to  compute 
a  new  Z  for  every  bed.  If  Z  is  the  value  of  Z  for  the 
inversion  horizon  the  value  of  Z(i+1)  for  the  overlying 
layers  is  given  by  equation  (10): 


' 
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Z(i  +  1)  =  Z  (i)  -  BTK(i-H)  (10) 

where  i  -  1,  2,  3...,  is  the  total  number  of  beds  overlying 
the  inversion  horizon,  and  BTK(i*1)  is  the  bed  thickness  of 
(i  +1)th  layer. 

After  the  gravitational  anomaly  of  all  layers  above  the 
inversion  horizon  are  calculated  and  summed  (  g  (sum)  of 
flowchart) ,  the  difference  between  the  observed  gravity 
(  g  (total)  of  flowchart)  and  the  gravity  g  (sum)  will  give  a 
new  version  of  the  gravity  g (inv)  which  can  be  used  for 
subsequent  iterations.  But,  if  the  initial  value  used  for 
the  first  iteration  is  too  low  the  inverted  h (x)  will  be  too 
small.  Since  g  (sum)  depends  heavily  on  h  (x) ,  keeping  other 
factors  unchanged,  too  small  a  h(x)  will  generate  too  low  a 
gravity  g  (sum)  .  Thus  the  difference  between  g  (total)  and 
g  (sum)  will  be  too  high,  and  the  inverted  h(x)  of 
subsequent  iteration  will  be  too  high  as  well.  Hence,  h(x) 
will  either  be  too  low  or  too  high  and  never  converge. 
Likewise,  if  the  initial  value  is  too  high,  h  (X)  will  never 
converge.  To  avoid  this  the  mean  of  g(res),  which  is 
(g  (total)  -  g(sum)  )  ,  and  the  gravity  g  (inv)  employed  for 

the  last  iteration  is  used  instead.  This  avoids  dependence 
on  a  judicious  choice  of  FRACTION,  but  also  avoids  the 
problem  of  divergence  resulting  from  a  wrong  choice  of 
FRACTION.  Since  the  gravity  used  for  every  iteration  is  not 
exactly  the  same,  at  least  for  the  first  few  iterations,  the 
rate  of  convergence  will  be  slower  than  that  using  a 


15 


Equo  II  on  (5) 


'  *4^ 


Figure  21  Schematic  representation  of  the 
multiioyer  inversion  procedure. 
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constant  gravity,  and  more  iterations  are  needed  to  reach 
convergence. 

2.3.  Oldenburg's  Cosine  Model 

To  test  the  one  layer  algorithm  of  Parker  and  Oldenburg 
it  was  decided  to  reproduce  some  of  the  results  obtained  by 
Oldenburg  (1974).  As  a  further  check,  Talwani's  (1959)  line- 
integral  method  for  two-dimensional  bodies  was  also  used  to 
compute  gravitational  attraction. 

The  numerical  model  used  by  Oldenburg  (1974)  was  a 
cosine  model,  4  km  high,  20  km  wide,  and  floored  7  km  deep 
(Figure  2.2A).  The  shape  of  this  model  is  given  by: 

h(x)  =2+2  Cos  (trX/10)  - 1 0<x<  1 0  (11) 

h (x)  =  0  j  x  |  > 10 

The  gravity  anomaly  (Figure  2.2C)  of  this  model  was 
computed  at  128  points,  each  1  km  apart,  using  Parker's 
algorithm.  The  parameters  used  were  Z  =  7 

km.  d  =1  gm/cc  and  H/Zo  =  4/7  which  is  0.57.  Talwain's  line- 
integral  method  was  also  used  to  calculate  the  gravity 
anomaly.  The  difference  between  the  two  resultant  profiles 
was  less  than  a  milligal. 

This  gravity  profile  (Figure  2.2C)  was  then  used  to 
find  the  shape  of  the  inversion  model.  Convergence  criteria 
were  the  same  as  that  used  by  Oldenburg.  The  inversion 
parameters  were  d=  1  gm/cc,  k(WH)  =  0.075,  k  (SH)  =  0.125  and 


Figure  2.2:  Two-dimensional  cosine  model  and  its 
gravity  anomaly.  Fig  2.2A  shows  the  normal  position  of 
the  origin  plane,  z  =  7  km,  and  the  position  for  optimum 
convergence,  z  —  5  km.  Fig  2.2B  and  Fig  2.2C  are  the 
gravity  profiles  computed  with  z  =  5  km  and  7  km  respectively. 
Note  that  the  vertical  axis  is  doubly  valued  to  indicate 
vertical  distances  relative  to  the  initial  and  shifted 
origins . 


DEPTH  IN  KM.  GRAVITY  IN  MGAL 


17 


18 


Z  was  varied  from  2  km  to  7  km  at  an  increment  of  1  km.  The 
frequencies  were  chosen  because  they  were  the  ones  used  by 
Oldenburg  (1974,  p5 33)  to  reproduce  the  topography. 

Inversion  with  the  above  stated  parameters  or  similar 
ones  did  not  produce  convergence  at  any  level  of  z  .  Later 
investigation  indicated  that  an  offset  of  -84  milligals  was 
actually  added  to  the  gravity  anomaly  in  Figure  2.2C  before 
inversion.  The  reason  for  the  offset  was  that  during  his 
computation  of  the  gravity  anomaly  of  the  cosine  model 
Oldenburg  actually  subtracted  2  km  from  each  of  the  128 
points  describing  the  elevation  of  the  model.  Such  an 
operation  merely  shifts  the  origin  plane  up  2  km  (Figure 
2.2A)  and  does  not  change  the  physical  relationship  between 
the  model  and  the  datum  plane  Z  .  A  displacement  of  the 
origin  plane  does  not  affect  the  validity  of  equation  (1) 
but  it  does  alter  the  numerical  values  of  Z  ,  h(x),  and  H  as 
well.  Thus  h  (x)  is  negative  for  5  <  | x |  <  10,  Z  is 
decreased  from  7  km  to  5  km,  and  H/Z  changes  from  4/7  to 
2/5,  because  H  =  raax|h(x) |  =  2  km 

To  shift  the  origin  plane  up  by  2  km  is  equivalent  to 
subtracting  the  gravity  effect  of  a  Bouguer  slab  of  2  km 
thick  from  the  initial  gravity  waveform  (Hinson,  1976. 
Appendix  C) . 

Each  of  the  128  gravity  stations  was  therefore  offset  - 
84  milligals,  and  the  resulting  gravity  profile  (Figure 
2.2B)  was  used  for  inversion.  The  inversion  parameters  used 
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were  d=  1 

gm/cc. 

k  (WH) 

=  0.075, 

k(SH)  =  0.125 

and 

z 

was 

varied  from 

2  km 

to 

7  km  at 

an  increment 

of 

1 

km . 

Convergence 

occurred  at 

all  levels 

of  Z  ,  except 

for 

z 

=  6 

km  and  7  km.  For  Z  =  6  km  and  7  km*  the  scheme  diverged 
after  three  iterations.  The  right  inversion  model  (Z  =5 
km),  shown  in  Figure  2.3,  took  four  iterations  to  converge. 

To  see  why  initial  inversion  with  the  gravity  of  Figure 
2.2C  diverged,  and  why  an  offset  gravity  (Figure  2. 2B) 
produced  convergence,  it  is  necessary  to  understand  how 
different  parameters  control  convergence.  Convergence  of 
Pakrer's  algorithm  is  governed  by  d  ,Z  and  H/Z.  For  the 
inversion  scheme  convergence  is  also  controlled  by  the 
filter  parameters. 

The  ratio  H/Z  only  affects  the  infinite  sum  of  the 
Par ker-Oldenbur g  algorithm.  As  long  as  H/Zo  is  less  than  1 
or  not  much  greater  than  1  convergence  of  the  infinite  sum 
is  assured.  In  our  present  situation  d  has  no  effect  upon 
convergence  because  the  value  used  for  inversion  was  the 
same  as  that  used  in  the  computation  of  the  gravity  profile. 
As  already  indicated,  for  any  given  gravity  profile  there 
exists  a  maximum  value  for  Z  ,Z  (max) ,  for  every  set  of 
inversion  parameters  being  used;  beyond  which  the  iterative 
procedure  no  longer  converges.  If  Z  (max)  is  high  enough  to 
enclose  the  appropriate  value  of  Z  needed  to  produce  a 
successful  inversion  model  there  is  no  problem.  However,  if 
Z  (max)  is  small  enough  to  exclude  the  proper  value  of  Z  the 
inversion  algorithm  will  never  converge  even  when  the  right 
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Z  is  used. 

It  appears  that  for  the  gravity  profile  in  Figure  2.2c, 
and  using  k(WH)  =  0.075  and  k(SH)  =  0.125,  the  value  of  Z  (max) 
was  2  km.  This  is  because  the  inversion  scheme  only  tended 
to  converge  at  Z  =  2  km.  With  k(SH)  reduced  from  0.125  to 
0.100  the  inversion  algorithm  not  only  converged  at  Z  =  2 
km,  but  also  tended  to  converge  at  Z  =  3  km.  Thus  Z  (max)  was 
then  3  km.  Since  the  right  value  of  Z  to  produce  the  proper 
inversion  model  is  7  km,  which  is  off  the  above  stated 
ranges  of  Z  (max) ,  the  initial  inversion  procedure  therefore 
diverged. 

It  can  be  observed  that  the  effect  of  increased 
filtering  is  twofold.  Firstly,  for  any  given  Z  increased 
filtering  decreases  the  number  of  iterations  required  to 
deliver  convergence.  Secondly,  increased  filtering  raises 
the  value  of  Z  (max)  .  In  theory  the  values  of  k(WH)  and  k(SH) 
can  be  decreased  to  such  an  extent  as  to  raise  the  magnitude 
of  Z  (max)  to  enclose  the  Z  necessary  to  produce  the  right 
inversion  model.  However,  in  practice,  severe  filtering  is 
most  undesirable.  Severe  filtering  cuts  off  high  frequency 
components  necessary  to  form  sharp  discontinuities,  and  the 
inversion  model  loses  most  of  its  detail.  The  gravity 
computed  from  the  inversion  model  might  not  agree  with  the 
observed  gravity. 

Instead  of  reducing  the  magnitudes  of  k(WH)  and 
k(SH)  to  increase  the  range  of  Z  (max)  it  was  found  that  the 
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same  effect  can  be  achieved  by  adding  an  appropriate  amount 
of  offset  to  the  initial  waveform.  This  was  evidenced  by  the 
inversion  using  the  gravity  profile  shown  in  Figure  2.2B. 
Using  this  gravity  waveform  the  iterative  scheme  converged 
at  all  levels  of  Z  except  at  Z  =  6  km  and  7  km.  Z  (max)  was 
extended  to  5  km. 

By  introducing  an  appropriate  amount  of  gravity  offset 
the  value  of  Z  needed  to  produce  the  right  inversion  model 
is  also  reduced.  With  the  cosine  model  at  a  base  of  7  km  (no 
offset  applied  to  the  initial  gravity)  Z  has  to  be  equal  to 
7  km  to  produce  the  right  shape.  However,  with  the  cosine 
model  based  at  2  km  below  the  origin  plane  Z  only  has  to  be 
equal  to  5  km  to  reproduce  the  input  shape.  As  already 
indicated  decreasing  Z  increases  the  rate  of  convergence. 
Hence,  fewer  iterations  are  needed. 

The  introduction  of  a  proper  amount  of  gravity 
displacement  also  has  the  effect  of  reducing  the  amount  of 
filtering  needed  to  produce  convergence.  To  decrease 
filtering  will  mean  allowing  high  wavenumber  oscillations  to 
pass.  This  is  especially  important  when  we  are  inverting 
bodies  with  sharp  discontinuities,  because  high  wavenumbers 
are  needed  to  form  sharp  curves.  To  achieve  fast  convergence 
of  the  inversion  scheme  it  is  most  desirable  if  Z  and  H/Z 
are  as  small  as  possible.  However,  decreasing  the  value  of 
Z  also  increases  the  ratio  H/Z  .  Hence,  the  infinite  sum 
takes  more  terms  to  converge.  If  Z  is  so  small  that  H/Z  is 
close  to  one  or  greater  than  one,  the  infinite  sum  will 
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DISTANCE  IN  KH 


OISTRNCE  IN  KH 


Fiqure  2.4:  Inversion  model  obtained  from  the  offset 
ravity  anomaly  (-126  ragal)  of  the  cosine  model.  Inversion 
arameters  used  were  4  km,  d-  1  qm/cc ,  ^  -0.075,  k 
.125.  iterations  were  used. 


SH 


Figure  2.3:  Inversion  model  obtained  from  the  offset 
gravity  anomaly  (-84  mqal)  of  the  cosine  model.  Inversion 
parameters  used  were  z  ■  5  km,  d-  1  qm/cc,  k^  =0.075,  kc  ■ 
0.125.  4  iterations  were  used.  SH 
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Fiouro  2.5:  Inversion  model  obtained  from  the  offset 
gravity  anomaly  (-126  maal)  of  the  cosine  model.  Inversion 
I aram» tors  wore  the  same  as  those  used  in  Fig  2.4  except 
K^O.125,  kc; f f  0.25,  4  iterations  were  used. 


Fiqure  2.6:  Difference  between  observed  gravity  ind 
gravity  calculated  from  the  inversion  models  shown  ir.  Fig 
2.3,  Fig  2.4,  and  2.5. 
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either  converge  very  slowly  or  diverge.  On  the  other  hand, 
if  Z  is  large  H/Z  will  be  small.  A  large  Z  will  mean  more 
iterations.  The  amount  of  computer  time  that  might  have  been 
saved  by  using  a  small  Z  could  well  be  cancelled  by  the 
increase  of  computer  time  in  calculating  the  extra  terms  of 
the  infinite  sum.  A  compromise  between  the  values  of  Z  and 
the  ratio  H/Z  is  therefore  needed  to  facilitate  convergence. 

Figure  2.4  is  the  inversion  model  produced  by  adding 
126  mgals  (3  km  Bouguer  Slab)  to  the  gravity  in  Figure  2.2C. 
With  Z  =  4  km  and  using  the  same  filter  (k(WH)=  0.075  and 

k  (SH) =  0.125)  it  took  three  iterations  to  deliver 

convergence.  When  the  filter  parameters  were  changed  to 
k(WH)  =  0.15,  k  (SH)  =  0.25  it  took  four  iterations  to 
converge.  Notice  the  shape  is  less  oscillatory  (Figure  2.5). 

Figure  2.6  shows  the  difference  between  the  observed 
gravity  and  the  gravity  obtained  from  the  inversion  models 
shown  in  Figure  2.3,  Figure  2.4  and  Figure  2.5. 

2.4.  Multilayer  Numerical  Example 

To  test  out  the  multilayer  inversion  program  the 
gravitational  anomaly  caused  by  a  four-layer  model  composed 
of  three  concentric  triangular  prisms  (Figure  2.7)  is  used. 
Each  triangular  prism  has  a  base  40  km  wide,  a  height  of 
4km,  and  a  density  contrast  of  0.5  gm/cm3.  The  lower 
triangular  model  is  based  at  8  km,  while  each  of  the 
overlying  layers  is  vertically  displaced  by  1  km. 
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The  gravity  of  this  four-layer  model  was  computed  for 
128  points,  each  1  km  apart,  using  the  equation  for 
triagular  prisms  in  Heiland  (p151,  1940,  1963).  The  gravity 
waveform  is  shown  in  Figure  2.8.  To  optimize  convergence 
251  milligals  were  added  to  each  of  the  128  gravity 
stations.  The  resultant  gravity  profile  was  then  used  for 
inversion.  The  parameters  used  for  inversion 

were:d  =  0.5  gm/cm3,  Z  =4  km.,  FRACTION  =  1/3, 

k(WH)  =  0.075,  k  (SH) =  0.125.  Six  iterations  were  used  to 
produce  convergence.  It  is  to  be  noted  that  by  adding  -251 
milligals  to  the  original  gravity  waveform  we  are  merely 
introducing  a  Bouguer  Slab  of  4  km  thick,  because  only  one- 
third  of  the  offset  gravity  is  used  for  inversion.  The 
inversion  model.  Figure  2.9,  agrees  well  with  the  input 
model.  Difference  between  the  observed  gravity  and  the 
inverted  gravity  is  shown  in  Figure  2.10.  Figure  2.11  shows 
the  difference  between  the  topography  of  the  input  model  and 
that  of  the  inversion  model. 

As  can  be  seen  in  Figure  2.10  (also  in  Figure  2.11)  the 
major  discrepencies  between  original  gravity  profile  (or 
topography)  and  inverted  gravity  waveform  (or  topography) 
occur  at  locations  with  sharp  discontinuity.  This  is  because 
the  high  frequency  components  which  are  needed  to  form  sharp 
curves  were  filtered. 
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Figure  2.7:  Four-layer  numerical  model,  composed  of 
three  concentric  triangular  prisms. 


Figure  2.8:  Original  gravity  and  gravity  computed  from 
the  inversion  model  of  the  four-layer  numerical  model. 


figure  2.9:  Comparison  of  the  inverted  shape  and  the 
■  riginal  shape  of  the  four-loyci  model. 


Fiqure  2.10:  Piffercnce  between  observed  gravity  \nd 
inverted  gravity  of  the  four- layer  model. 


DIFFERENCE  BETWEEN  ORG.  AND  INV.  TOPOG. 
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Figure  2.11:  Difference  between  the  inverted  topography 
and  the  original  topography  of  the  four-layer  model. 
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CHAPTER  3  GEOLOGY  AND  GEOPHYSICS  OF  THE  AREA 


3.1.  Geological  Background 

The  locations  of  the  two  gravity  profiles  for  which  we 
are  going  to  do  the  inversion  are  shown  in  Figure  3.1.  This 
area  is  of  great  economical  and  geological  interest. 
Significant  amounts  of  oil  and  gas  have  accumulated  in  the 
Phanerozoic  sediments  that  overlie  the  basement.  Nickel  is 
being  mined  from  the  Thompson  nickel  belt,  one  of  the 
world's  major  nickel  deposits.  The  area  is  also  of  tectonic 
interest  because  the  boundary  zone  (Bell,  1971b)  between 
Churchill  and  Superior  provinces  which  outcrops  in  the 
northern  part  of  Manitoba  might  well  extend  south  (Bell, 
1971b)  beneath  the  Phanerozoic  sediments  (Figure  3.1). 

The  southern  plains  of  Manitoba  and  Saskatchewan  form 
part  of  the  Interior  Platform.  Phanerozoic  sedimentary  rocks 
thin  out  towards  eastern  Manitoba.  The  nearly  horizontal 
sedimentary  rocks  are  covered  by  a  thick  mantle  of  glacial 
drift  to  form  the  plains  and  plateaux  of  the  Interior  Plains 
physiographical  province.  Beneath  the  Phanerozoic  strata  are 
Precambrian  sedimentary  and  crystalline  rocks,  the  latter 
forming  the  basement  of  this  area.  To  the  north  of  the 
plains  the  Precambrian  rocks  are  exposed  on  the  Canadian 
Shield.  The  Precambrian  rocks  of  northern  Saskatchewan  and 
northwest  Manitoba  form  part  of  the  Churchill  structural 
province,  while  those  of  north-east  Manitoba  compose  a 
portion  of  the  Superior  structural  province.  Part  of  the 
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Figure  3.1:  Gravity  map  of  the  southern  plains  of 
Saskatchewan  and  Manitoba,  showing  locations  of  the  two 
gravity  profiles,  and  area  of  detail  gravity  survey. 
SU=Suf field,  SC=Swif t  Current,  MH=Medicine  Hat,  E=Elbow, 
R= Regina,  and  W=Winnipeg. 
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Superior  Province  is  also  exposed  in  southeast  Manitoba,  The 
Superior  and  Churchill  provinces  are  partially  separated  by 
a  transition  zone:  the  Pikwitonei  province  (Bell,  1971b). 

3.2.  Regional  Geology  and  Geophysics  of  the  Area 

Superior  Province 

The  Manitoba  portion  of  the  Superior  province  is 
dominated  by  east-west  trending  structures.  Its  rocks 
consists  of  belts  of  calc-alkaline  and  calcic  volcanic 
sequences  with  interlayered  sediments,  and  belts  of 
sedimentary  gneiss  interlaminated  with  granitic  intrusions 
(Wilson,  1971).  The  volcanic-sedimentary  belts  have  been 
extensively  intruded  by  granitic  plutons,  and  have  been 
deformed  and  metamorphosed  to  the  greenschist  metamorphic 
facies.  The  sedimentary  gneisses  have  been  metamorphosed  to 
amphibolite  and  occasionally  granulite  facies. 

Strong  east-west  magnetic  and  gravitational  trends  are 
observed.  The  two  types  of  belt  are  characterised  by 
distinct  magnetic  patterns.  On  the  whole  it  appears  as  belts 
of  magnetically  high  areas  composed  of  broad  linear 
anomalies  of  moderate  relief,  and  narrow  linear  anomalies  of 
high  relief.  The  volcanic-sedimentary  belts  occur  in  areas 
of  lower  magnetic  intensity  than  the  surrounding  gneissic 
areas  which  are  broad  E-W  trending  magnetic  highs  with 
little  relief.  Local  sharp  and  intense  magnetic  highs  are 
registered  across  the  volcanics  and  ironstone  formations. 
Granitic  bodies  are  characterised  by  elliptical  high 
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amplitude  anomalies  (Kornick,  1971;  Kornick  and  Maclaren, 
1966;  Hall ,  1968,  1971,  1976). 

The  Boundary  Zone 

There  is  much  controversy  concerning  the  exact  location 
and  the  extent  of  the  boundary  between  the  Churchill  and 
Superior  structural  provinces.  For  the  past  several  decades, 
earth  scientists  have  tried  to  use  both  geological  and 
geophysical  data  to  infer  the  exact  location  and  extent  of 
this  boundary,  however,  the  boundary  enigma  is  still 
unsolved.  A  transition  zone,  rather  than  a  sharp  boundary, 
was  proposed  to  separate  the  Churchill  and  Superior 
provinces  (Burwash  et  al.  1962;  Peterman  and  Hedge,  1964; 
Bell,  1971b) .  Geochronological  data  from  bottom-hole  samples 
indicated  that  the  Churchill  and  Superior  provinces  of  the 
Canadian  Shield  not  only  extend  into  southern  Saskatchewan 
and  Manitoba  but  also  into  the  Williston  basin  of  North 
Dakota  (Burwash  et  al.  1962;  Peterman,  1962;  Peterman  and 
Hedge,  1964) .  On  the  basis  of  gravity  data  (Innes,  1960; 
Wilson  and  Brisbin,  1961)  and  K-Ar  ages  of  subsurface 
samples  (Burwash  et  al.  1962;  Peterman  and  Hedge,  1964),  it 
was  indicated  that  the  boundary  zone  in  the  subsurface 
extends  from  the  Nelson  River  gneissic  zone  southwest  into 
eastern  Saskatchewan  and  hence  into  north-central  North 
Dakota,  following  the  Nelson  River  gravity  low.  North  of 
Lake  Winnipeg,  Bell  (1971b)  defined  the  boundary  zone  as  the 
Pikwitonei  province,  which  can  be  traced  intermittently  from 
the  foot  of  Lake  Winnipeg,  northeast,  then  east,  under  the 


. 
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Hudson  Bay  Lowland  and  the  flat-lying  Aphebian  rocks  that 
outcrop  west  of  James  Bay  (Bell,  1971b,  p32.  Figure  8).  The 
rocks  north  and  northwest  of  the  Pikwitonei  province  -  the 
Waboden  subprovince,  were  initially  included  by  Bell  as  part 
of  the  Churchill  province.  However,  recent  geological  field 
work  and  Rb-Sr  radiometric  dating  by  Cranstone  and  Turek 
(1976)  have  shown  that  the  Wabowden  subprovince  is  best 
described  either  as  a  separate  entity  or  as  part  of  the 
Pikwitonei  province. 

The  Pikwitonei  province  is  a  zone  of  early  Precambrian, 
retrograded,  granulite  to  amphibolite  facies  gneisses  that 
lie  uncomf ormably  beneath  'type*  Archean  greenstone-belt 
rocks  -  Cross  Lake  Group  (Bell  1971).  It‘s  rocks  have  a  NE- 
SW  fabric  which  abruptly  truncated  the  E-W  trends  of  the 
Superior  Province. 

Like  the  Pikwitonei  province,  the  Waboden  subprovince 
is  NE-SW  trending.  Its  rock  types  are  predominantly  migmatic 
gneisses  which  enclose  a  belt  of  sedimentary  and  volcanic 
rocks  of  significantly  lower  metamor phic  grade.  In  general, 
the  metasedimentary  and  metavolcanic  rocks  appear  to  be 
greenschist  to  lower  amphibolite  facies,  whereas  the  layered 
migmatic  gneisses  are  generally  upper  amphibolite  facies 
(Bell,  1971b).  The  vast  nickel  deposits  of  the  ‘Thompson 
Nickel  Belt*  occur  in  the  west  edge  of  this  subprovince.  The 
boundary  between  the  Wabowden  subprovince  and  the  Pikwitonei 
province  is  a  major  fault  zone  (Assean  Lake  fault  zone)  in 
the  north  (Bell,  1971b)  and  a  metamorphic  transition  zone  in 
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the  south  and  central  regions  (Cranstone  and  Turnek,  1976). 
The  definition  of  the  boundary  in  the  south  is  not  well 
established  as  similar  gneisses  occur  in  both  provinces. 

Aeromagnetic  and  gravity  studies  indicated  that  there 
is  a  remarkable  correlation  of  surface  geology  with  gravity 
and  magnetic  anomalies  (Bell,  1971b)  in  the  Pikwitonei 
province.  The  Pikwitonei  province  granulite  facies  gneisses 
are  reflected  as  small  ovoid  high  and  low  magnetic  anomalies 
of  relatively  high  amplitude  (Bell,  1971b;  Kornik  and 
Maclaren,  1966;  Kornik,  1969).  In  contrast  with  the 
Pikwitonei  gneisses  the  Waboden  migmatic  gneisses  is 
characterised  by  elongated,  high  intensity  magnetic 
anomalies;  parallel  to  the  trend  of  the  belt.  The  Pikwitonei 
province  and  the  Waboden  subprovince  are  marked  by  a 
prominent  gravity  high  and  low  respectively,  of  similar 
magnitude  and  dimensions  (Wilson  and  Brisbin,  1961,  1962; 
Gibb  1968b) .  It  was  on  the  basis  of  these  regional 
geophysical  trends  that  Bell  extrapolated  the  boundary  zone 
into  the  areas  covered  by  Paleozoic  and  Pleistocene  deposits 
(Bell,  1971b,  p32.  Figure  8). 

The  eastern  limit  of  the  Churchill  province  proposed  by 
Bell  (1971b)  is  on  the  east  flank  of  the  Nelson  River 
gravity  high,  while  that  of  Innes(1960),  Wilson  and  Brisbin 
(1961),  Burwash  (1962)  and  Peterman  (1964)  is  along  the 
Nelson  River  gravity  low. 


Churchill  Province 
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Churchill  province  rocks  consist  of  mainly 
metasediments  and  migraatites,  plus  minor  metamorphosed 
intermediate  to  mafic  volcanic  rocks.  Bodies  of  massive  to 
weakly  foliated  granitic  rocks  are  also  present.  The 
volcanic-sedimentary  areas  in  the  Churchill  province  are 
different  magnetically  from  the  volcanic-sedimentary  areas 
in  the  Superior  province.  The  Churchill  volcanic-sedimentary 
areas  do  not  contain  the  large  gnessic  belts  which  produce 
the  large  magnetic  highs  present  in  the  Superior  province. 
Instead  it  is  dominated  by  local  intrusive  features  which 
have  smaller  and  discontinuous  magnetic  anomalies.  Like  the 
magnetic  anomalies,  the  gravity  anomalies  are  isolated  and 
closed  and  do  not  show  any  strong  trends. 

3.3.  Seismic  Data 

Much  seismic  work  have  been  done  to  determine  the 
crustal  structure  of  Western  Canada.  Seismic  surveys  carried 
out  by  G.L.  Cumming  and  E.R.  Kanasewich  et  al  (1966)  in  the 
plains  of  southern  Alberta  and  Saskatchewan  indicated  a 
crustal  thickness  of  about  47  km  in  eastern  Alberta,  and  43 
km  in  western  Saskatchewan.  Their  derived  crustal  model  from 
Swift  Current,  Saskatchewan  to  Suffield,  Alberta  is 
summarized  in  Table  I.  Similar  surveys  were  done  in  the 
southern  part  of  eastern  Manitoba  (Hall  and  Hajnal,  1969; 
Gurbuz,  1969,  1970).  Their  derived  crustal  thickness  is 
about  34  km.  The  crustal  model  adopted  by  Gurbuz,  (1969, 
1970)  is  shown  in  Table  II.  It  is  based  on  the  data  present 
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in  these  two  tables,  together  with  the  seismic  data 
(Hall, 1969)  concerning  the  structure  of  the  crust  around 
central-south  Lake  Winnipeg  that  the  parameters  (e.g.  depth 
of  interface,  density  and  thickness)  necessary  for  the 
multilayer  inversion  scheme  are  derived. 


' 


Table  I 


1 - - - 

■j - 

- - r 

-~i 

| Discontinuity |  SUFFIELD  | 

SWIFT  CURRENT 

1 

i   . 

« 

.  .  j 

i 

■■I 

j  Depth 

P-Wa ve 

Density 

Depth 

1 

1  (Km) 

Velocity 

(Km) 

1 

(Km/sec) 

(Gm/cc) 

1 

i 

|  Surface 

|  0.0 

0.0 

i 

1 

i 

i 

1 

|  Basement 

j  2.24 

2.02 

1 

1 

i 

I 

1 

6.16 

2.73 

I 

1 

|  Man.  Riel 

J  12.5 

10.4 

1 

i 

i 

6.5 

2.8 

i 

i 

i 

|  Alta.  Riel 

|  35.5 

33.4 

i 

i 

7.15 

3.17 

i 

i 

i 

|  Mohorovicic 

S  46.9 

43.2 

i 

i 

1 

i 

8.08 

3.4 

i 

! 

i 

i 

d 

I 

.i 

Crustal  model  from 

Suf field. 

Alberta 

to  Swift 

Current,  Saskatchewan.  All  depths  are  measured 
relative  to  the  surface.  The  Table  is  not  drawn  to 
scale. 


Density  values  are  calculated  from  the 
density-compressiona 1  velocity  curve  given  by  Nafe 
and  Drake  (1957)  . 
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Table  II 
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Crustal  model  of  the  southern  portion  of 
eastern  Manitoba,  as  derived  by  Gurbuz  (1969, 
p970) .  All  depths  are  measured  relative  to  the 
surface.  All  density  values  are  computed  from  the 
density-compressiona 1  velocity  curve  given  by 
Nafe  and  Drake  (1957). 


37 


It  is  to  be  noted  that  there  is  ambiguity  in  the  usage 
of  the  term  "Riel  discontinuity"  (R.D.).  Previously 
Kanasewich  and  Curaming  et  al  (1959)  termed  the  discontinuity 
just  above  the  Mohovicic  discontinuity  (Moho)  the  "Conrad 
Discontinuity".  This  nomenclature  was  criticized  because  it 
implied  a  stratigraphic  or  petrologic  correlation  with  the 
European  Conrad  discontinuity  which  represents  an  interface 
with  velocity  changing  from  about  5.6  to  6.2  km/sec. 
(Conrad,  1925;  Jeffreys,  1926) .  In  southern  Alberta,  the 
refracting  and  reflecting  intermediate  layer  represents  a 
change  in  velocity  from  6.5  to  7.2  km/sec.  To  avoid  the 
implication  of  an  intercontinental  correlation,  Kanasewich 
and  Cumming  et  al  (1968)  subsequently  replaced  the  term 
"Conrad"  with  "Riel".  In  Manitoba  the  term  "Riel 
discontinuity"  was  used  to  refer  to  the  discontinuity  which 
represents  a  velocity  change  from  6.1  to  6.8  km/sec. 

To  avoid  any  confusion  in  the  future  we  designate  the 
Alberta  R.D.  (ARD)  as  the  discontinuity  just  above  the 
intermediate  layer  with  velocity  7.1  ±0.15  km/sec;  and 
Manitoba  R.D.  (MRD)  as  the  interface  just  above  the  layer 
with  velocity  6.72  ±  0.23  km/sec. 
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CHAPTER  4 


4.1.  Source  of  Data 

The  gravity  data  used  for  inversion  were  obtained  from 
the  Dominion  Observatory.  Additional  data  were  obtained  by 
C.  McCloughan  and  the  author  during  a  gravity  field  program 
in  the  summer  1976.  The  gravity  values  of  stations  marked 
along  the  two  profiles  were  used  as  the  gravity 
observations. 

The  Gravity  Survey  (Figure  4.1) 

In  the  summer  of  1976,  Mr.  McCloughan  and  the  author, 
under  the  direct  supervision  of  Professor  G.L.  Cumming  of 
the  University  of  Alberta,  undertook  a  gravity  survey  in  the 
southern  plains  of  Saskatchewan  and  Manitoba.  The  survey  was 
intended  to  add  more  details  to  the  Dominion  Observatory 
Gravity  Map  Series.  Observations  were  carried  out  with  a 
Lacoste  and  Romberg  gravimeter,  supplied  by  the  Dept.  of 
Energy,  Mines  and  Resources  (E.M.R.).  The  survey  area  is 
bounded  by  50°30'N  and  52°N  latitudes  and  101°W  and  104°W 
longtitudes.  About  690  observations  were  made.  Readings  were 
taken  at  1  or  2  miles  intervals,  depending  on  the  density  of 
gravity  stations  already  present  in  the  area.  Observations 
were  taken  along  highways,  railway  stations  and  crossings 
and  grid  roads.  Elevation  control  was  provided  by  known 
elevations  at  section  corners,  railway  crossings, 
benchmarks,  lake  levels  and  by  the  use  of  two  altimeters. 
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Readings  were  taken  at  gravity  control  stations  every  2 
hours  and  at  points  of  known  elevation  every  hour. 

The  station  positions  were  marked  on  4-mile  topographic 
mapsr  from  which  their  latitudes  and  longitudes  were  scaled. 
The  updated  version  of  the  Potsdam  Standard  gravity  base  was 
used  in  data  reduction.  The  free  air  and  Bouguer  correction 
amounts  to  0.1968  milligals  per  meter  above  sea  level.  A 
density  of  2.67  gms.  per  c.c.  was  used  in  correcting  for  the 
attraction  of  material  between  the  station  and  sea  level. 
Topographic  correction  was  not  applied  to  the  data  because 
there  is  not  much  variation  in  terrain  in  the  area.  The 
elevation  of  most  of  the  stations  obtained  are  believed 
accurate  within  5  feet.  The  estimated  uncertainty  of  the 
final  Bouguer  anomalies  is  ±0.3  mgal. 

It  is  to  be  noted  that  the  gravity  values  appeared  in 
the  above  stated  gravity  map  series  were  obtained  using  the 
old  version  of  the  Potsdam  Standard  for  gravity.  An  analysis 
of  the  Dominion  gravity  data  and  our  data  at  the  same 
locations  indicated  that  our  gravity  values  are  7  mgals 
lower  than  that  of  the  Dominion  Observatory.  To  bring  all 
Dominion  data  to  the  new  Potsdam  Standard  we  subtract  7 
mgals  from  all  Dominion  gravity  stations. 

Location  (Figure  3.1) 

Profile  (1)  runs  from  Suffield  i111o00*W)  Alberta, 
along  50°20,N  latitude,  to  the  northern  part  of  Winnipeg 
(96°00 • W) .  Profile  (2)  is  50.5  minutes  north  of  Profile  (1) 
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Figure  4.1:  Detail  gravity  map  of  southern 
Saskatchewan . 


4  1 


and  runs  along  latitude  51°10.5,N. 

f2£  the  Sediment  La  yer 

The  sedimentary  basin  of  the  southern  plains  of 
Saskatchewan  and  Manitoba  can  be  roughly  divided  into  three 
groups.  Sediments  with  a  mean  density  of  2.42  gm/cc  are 
considered  to  be  the  first  group,  and  they  are  composed  of 
shale,  shaly  sand,  sandstone  and  some  limestone.  The  second 
group  consists  of  evaporites  which  have  a  mean  density  of 
2.16  gm/cc.  Dolomites,  limestones,  high  velocity  sands  and 
shales  form  the  third  group,  and  they  are  assumed  to  have  a 
mean  density  of  2.67  gm/cc.  The  thickness  and  altitude  of 
various  formations  were  determined  from  isopach  maps1.  The 
resultant  gravitational  contribution  of  the  sediments,  which 
were  computed  using  Talwani*s  (1959)  method,  was 
subsequently  subtracted  from  their  corresponding 
gravitational  observations  on  the  profile  to  give  the 
reduced  gravitational  profile  for  inversion. 

4.2.  Crustal  Models 

In  the  following  inversions  four  different  crustal 
models  are  assumed,  and  they  are  shown  in  Figure  4.2. 

Model  I  (Fig .  4.2-A) 

Model  I  consists  of  a  sedimentary  layer  underlain  by 
three  crustal  layers.  A,  B  and  C  with  a  discrete 


i  From  Geological  History  of  Western  Canada,  ASPG,  Calgary, 
Alberta . 


' 
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compressiona 1  velocity. 

Model  II  (Fig.  4.2-B) 

In  this  model,  ARD  is  assumed  non-existent  and  the 
material  between  MRD  and  the  Moho  is  represented  by  a  layer 
(BC)  whose  density  is  the  mean  of  the  densities  of  layers  B 
and  C  in  Model  I.  The  expression  used  to  find  the  mean 

density  is: 

i=n  i=n 

Vd(i)h(i)/  £h(i)  n  =  1,2,3... 

fi  i=i  ~  total  number  of  layers 

where  d (i)  and  h(i)  is  the  density  and  thickness  of  the  ith 

layer  respectively. 

Model  III  (Fig.  4.2-C) 

This  model  is  similar  to  Model  II  except  MRD  is  assumed 
non-existent  instead  of  ARD.  The  mean  of  the  densities  of 
layers  A  and  B  of  Model  I  is  used  as  the  density  of  the 
layer  (AB)  between  ARD  and  the  basement. 

Model  IV  (Fig.  4.2-D) 

In  this  model,  layers  A,  B  and  C  of  Model  I  is  combined 
to  form  one  single  layer  (ABC)  whose  density  is  the  mean  of 
the  densities  of  layers  A,  B  and  C. 

4.3.  Inversion  of  Profile  (1) 

Profile  (1)  is  divided  into  two  segments:  east  of  Swift 
Current,  and  west  of  Swift  Current.  To  see  how  the  various 
inversion  models  produced  by  the  multilayer  inversion  scheme 
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compared  with  the  Suffield  -  Swift  Current  crustal  model  put 
forward  by  Cumming  and  Kanasewich  et  al  (1966) ,  inversion 
was  first  performed  on  the  western  segment. 

Suffield  -  Swift  Current 

We  used  as  gravitational  observations  the  first  35 
stations  east  of  111°58*N  longitude  iFig.  4.3).  Sediment 
gravity  at  corresponding  locations  was  computed  and  removed 
to  obtain  the  reduced  gravitation  profile.  In  order  to  meet 
the  requirement  that  for  finite  transform  the  data  series 
has  to  be  infinitely  replicated,  a  "mirror  image"  of  the 
reduced  gravitational  profile  was  generated.  A  spline 
routine  was  used  to  interpolate  these  stations.  The 
resulting  profile  was  then  digitised  at  intervals  of  5  km 
yielding  a  series  128  points  long  (Fig.  4.4).  Convergence 
criteria  were  the  same  as  the  numerical  examples.  A  maximum 
of  15  iterations  were  allowed. 

Table  III  shows  the  thicknesses,  densities  and  depths 
of  each  layer  at  Swift  Current  and  Suffield  (Cumming  and 
Kanasewich  1966).  It  is  to  be  noted  that  these  figures  are 
the  average  values  around  Swift  Current  and  Suffield. 
However,  during  the  process  of  inversion,  these  values  were 
used  as  the  values  of  the  corresponding  layer  at  Swift 
Current  and  Suffield.  It  was  assumed  that  the  thickness  (or 
density)  of  a  layer  vary  linearly  along  the  profile  if  the 
thickness  (or  density)  of  the  stratum  is  not  the  same  at 
both  ends  of  the  profile. 
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The  total  number  of  iterations  required  for 
convergence,  and  the  inversion  parameters  used  for  the  four 
different  models  are  also  shown  in  Table  III.  The  depth  of 
the  Moho  at  Swift  Current  provides  the  constraints  required 
for  the  proper  choice  of  Z  •  It  was  assumed  that  the  depths 
of  the  Moho,  ARD,  and  MRD  at  Swift  Current  are  at  43.2,  33.4 
and  10.4  km  respectively. 

The  above  described  procedure  and  assumptions  also 
apply  to  the  inversions  that  follow. 


Table  III 


MODEL  I 


|  INVERSION 

|  PARAMETERS 

I - — - 

|  LOCATION 

| - - r 

|  Average  I 

(  Thickness  of  | 

(  Layer  (km)  | 

I - - - — 4 

j  Average  | 

|  Density  I 

I  of  | 

|  Layer  (gm/c.c.)  j 

I - L 


Z  =33.5  km 
k  =0.087 
k  =0.057 
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- ( - - 

A  |  10.26 

B  I  23.00 
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B 

C 


2.7  3 
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z 

k 
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A 

BC 
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MODEL  III 
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Table  III.  Suf field  -  Swift  Current  inversion  parameters. 
ITER  =  number  of  iterations  required  for 
convergence. 
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Figure  4.3:  Gravitational  observations  of  the  Swift 
Current-Suf field  profile. 
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Figure  4.4:  Reduced  gravitational  observations  of  the 
Swift  Current-Suf field  profile. 
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Figure  4.5:  Inversion  model  of  Swift  Current-Suf field 
profile  assuming  crustal  model  I.  Vertical  exaggeration 
=  7. 
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Figure  4.6:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.5. 
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Figure  4.7:  Inversion  model  of  Swift  Current-Suf field 
profile  assuming  model  II. 
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INVERTED  AND  ORIGINAL 
GRAVITY  ANOMALIES 


Figure  4.8:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.7. 
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INVERTED  TOPOGRAPHY 


Figure  4.9:  Inversion  model  of  Swift  Current  Suf field 
profile  assuming  crustal  model  III. 
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INVERTED  AND  ORIGINAL 
GRAVITY  ANOMALIES 


Figure  4.10:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.9. 
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Figure  4.11:  Inversion  model  of  Swift  Current-Suf f ield 
profile  assuming  crustal  model  IV. 
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Figure  4.12:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.11. 
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The  inversion  models  and  the  gravity  profiles  from 
inversion  for  the  four  different  crustal  models  are  shown  in 
Fig.  4.5,  4.6,  4.7,  4.8,  4.9,  4.10,  4.11  and  4.12. 

Mo del  I  (Fig.  4.5  and  Fig.  4.6) 

The  inverted  topography  of  this  model  appears  to  be  in 
reasonable  agreement  with  that  of  Gumming  and  Kanasewich  et 
al  (1966).  The  average  inverted  depth  of  the  Moho  around 
Suffield  (from  45  to  60  km  in  the  figure)  is  45.7  km  , 
compared  with  46.9  km  from  seismic  data.  The  average 
inverted  depths  of  the  ARD  (34.5  km)  and  MRD  (11.5  km)  also 
agree  reasonably  well  with  those  obtained  by  Cumming  and 
Kanasewich  et  al  (1966,  Table  III). 

Seismic  determinations  (Cumming  and  Kanasewich  1966) 
indicated  that  east  of  Suffield  the  crust  thins  out 
gradually  towards  Swift  Current.  Such  trend  can  also  be  seen 
in  the  inversion  model.  Though  there  is  no  clear  seismic 
indication  of  the  relief  of  the  Moho  between  70  and  140  km, 
this  does  not  mean  the  inverted  topography  of  this  area  is 
unreal.  It  is  just  that  there  was  insufficient  seismic  data 
available. 

The  agreement  between  the  observed  and  the  anomaly 
calculated  from  the  inversion  model  is  shown  in  Figure  4.6. 
Models  II,  III  and  IV  (Fig.  4.7,4.8,4.9,4.10,4.11  and  4.12) 

Generally  speaking,  the  inverted  topography  of  the  Moho 
of  these  models  is  of  similar  shape  as  that  of  Model  I;  and 
all  3  models  show  crustal  thinning  towards  Swift  Current. 


< 
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The  inverted  depths  of  the  Moho  and  MRD  around  Suffield  of 
Model  II  appear  to  be  in  reasonable  agreement  with  that 
obtained  from  seismic  determination  (Cumming  and  Kanasewich, 
1966).  The  same  is  true  for  the  inverted  depths  of  the  Moho 
and  ARD  of  Model  III  and  the  Moho  of  Model  IV. 

A  comparison  of  the  inverted  shape  of  the  Moho  of  the  4 
models  shows  that  Models  I  and  II,  which  exhibit  similar 
relief,  have  the  lowest  relief  of  all  models;  followed  by 
Models  III  and  IV  in  ascending  order.  The  reason  why  Model 
IV  and  III  have  relatively  higher  relief  is  obvious:  to 
account  for  a  given  gravity  profile  the  topographic  relief 
of  a  structure  at  greater  depth  has  to  be  higher  than  that 
of  a  shallow  structure.  The  agreement  between  the  observed 
gravity  and  the  anomalies  calculated  from  the  3  inversion 
models  is  shown  in  Fig.  4.8,  4.10  and  4.12. 

On  the  whole,  the  anomalies  computed  from  all  4 
inversion  models  fit  the  original  observations  quite  well. 
The  short  wavelength  anomalies  probably  arise  from  error  in 
using  isopach  maps  to  compute  the  sediment  gravity;  the  non- 
two-dimensionality  of  the  sediment^cr ustal  boundary  and  the 
effect  of  possible  additional  shallow  geological  features, 
may  also  contribute  significantly  to  these  anomalies. 

Swift  Current  -  Winnipeg 

The  observed  gravity  of  this  segment  of  Profile  (1)  is 
shown  in  Figure  4.13;  while  the  reduced  gravitational 
profile  fs  shown  in  Figure  4.  14.  Procedure  similar  to  that 
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used  on  the  western  segment  was  employed  here  as  well.  The 
data  series  had  256  points;  at  a  sampling  interval  of  7  km 
Inversion  parameters,  together  with  densities,  depths  and 
thicknesses  of  each  layer,  used  for  the  four  different 
models  are  shown  in  Table  IV. 

It  can  be  observed  from  Table  IV  that  the  density  of 
layer  B,  and  the  thicknesses  of  various  layers  at  Swift 
Current  are  different  from  those  at  Winnipeg.  Linear 
interpolation  was  employed  to  find  the  densities  and 
thicknesses  of  intermediate  stations.  Again  the  depth  of  the 
Moho  at  Swift  Current  provides  the  necessary  contraint 
required  for  choosing  the  right  Z  .  The  depths  of  the  Moho 
ARD,  and  MRD  at  Swift  Current  are  assumed  at  43.2,  33.4  and 
10.4  km  respectively. 
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Table  IV.  Swift  Current 
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Figure  4.13:  Gravitational  observation 
Current-Winnipeg  profile. 


of  the  Swift 


Winnipeg 


GRAVITY  IN  MGAL 


65 


-100+ 

0 


T . r — — r  'r 

160  360  540  720 

DISTANCE  IN  KN 


900 


Figure  4.14:  Reduced  gravitational  observation 
the  Swift  Current-Winnipeg  profile. 
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Figure  4.15:  Inversion  model  of  Swift  Current-Winnipeg 
profile  assuming  crustal  model  I.  Vertical  exaggeration  =18 
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Figure  4.16:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.15. 
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Figure  4.17:  Inversion  model  of  Swift  Current-Winnipeg 
profile  assuming  crustal  model  II. 
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Figure  4.18:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.17. 
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Figure  4.19:  Inversion  model  of  Swift  Current-Winnipeg 
profile  assuming  crustal  model  III. 
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Figure  4.20:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.19. 
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Figure  4.21:  Inversion  model  of  Swift  Current-Winnipeg 
profile  assuming  crustal  model  IV. 
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Figure  4.22:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.21. 
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The  inversion  models  and  their  gravitational  profiles 
are  shown  in  Figures  4.15,  4.16,  4.17,  4.18,  4.19,  4.20, 
4.21  and  4.22. 

Model  I  (Fig.  4.15  and  Fig.  4.16) 

This  inversion  model  exhibits  a  gradual  decrease  in 
crustal  thickness  towards  Winnipeg.  The  mean  rate  of 
thinning  is  (0.013).  The  inverted  depths  of  the  MR D  and  ARD 
and  Moho  around  Winnipeg  are  20.5,  28  and  36.5  km 
respectively;  compared  with  the  average  depth  of  18±3  km  for 
the  MRD,  25.5±3.5  km  for  the  ARD  and  34.0±3  km  for  the  Moho 
(Gurbuz  1969,  1970).  The  Moho  shows  considerable  relief  from 
260  km  to  540  km.  The  ridge  around  470  km  corresponds  with 
the  southern  extension  of  the  "Nelson  River  Bouguer  gravity 
high".  The  east  flank  of  this  ridge  corresponds  with  the 
Churchill  -  Superior  boundary  proposed  by  Bell  (1971). 

Figure  4.16  shows  the  gravity  calculated  from  inversion 
fits  the  original  gravitational  observations  rather  well. 

Model  II  (Fig.  4.17  and  Fig.  4.18) 

This  model  exhibits  similar  topographic  relief  as  that 
obtained  in  Model  I.  The  crustal-thinning  trend  which  was 
indicated  in  Model  I  is  also  revealed  in  this  model. 
However,  the  mean  rate  of  crustal-thinning,  which  is  0.015, 
is  greater  in  this  model.  The  inverted  depths  of  the  MRD  and 
Moho  around  Winnipeg  are  18  and  34  km  respectively.  These 
values  agree  well  with  that  obtained  by  seismic 
determination  (Gurbuz  1967,  1970). 
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The  agreement  between  the  observed  gravity  and  the 
anomalies  calculated  from  the  inversion  model  is  shown  in 
Figure  4,18. 

Model  III  and  IV  (Fig.  4.19  to  4.22) 

Even  though  the  inverted  topography  of  the  Moho  of 
these  two  models  is  of  similar  shape  as  that  obtained  in 
Model  I  and  II,  however,  there  are  significant  differences. 
Firstly,  the  structural  relief  is  considerably  stronger  in 
these  two  models,  a  phenomenon  which  is  consistent  with  that 
observed  in  Model  III  and  IV  of  the  Suffield  -  Swift  Current 
profile.  Secondly,  though  these  two  models  still  exhibit 
crustal-thinning  towards  Winnipeg,  the  rate  of  thinning  is 
much  lower  than  that  in  Model  I  and  II.  Thirdly,  the 
inverted  depths  of  the  ARD  of  Model  III,  and  the  Moho  of 
Model  III  and  IV  around  Winnipeg  fail  to  reach  the  average 
depths  of  these  two  interfaces  determined  by  seismic  data 
(Gurbuz  1967,  1970) .  Around  Winnipeg  the  inverted  depths  of 
the  ARD  and  Moho  for  model  III  are  32  km  and  43  km 
respectively;  compared  with  25.5±3.5  km  for  the  ARD  and 
34 .0±3. 0  km  for  the  Moho  (Gurbuz  1967,  1970).  For  Model  IV 
the  inverted  depth  of  the  Moho  is  44  km  which  is  about  10  km 
deeper  than  the  average  value  determined  by  Gurbuz  (1967, 
1970)  . 

Although  the  two  gravitational  waveforms  (Fig  4.20  and 
4.22)  calculated  from  these  two  models  fit  the  observed 
gravitational  anomalies  well,  we  do  not  consider  these  two 
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models  acceptable  because  the  inverted  depths  of  the  ARD  and 
Moho  around  Winnipeg  do  not  agree  with  that  determined  by 
seismic  surveying  (Gurbuz  1967,  1970). 

4.4.  Inversion  of  Profile  (2) 

This  profile  runs  along  latitude  51°10.5'N  and  extends 
from  longitude  107°50*W  to  96°00*W.  It  will  be  referred  to 
as  the  "Elbow  (longtitude  106°35'W,  latitude  51°13*N)  *  Lake 
Winnipeg”  profile.  Seismic  data  (Gurbuz  1969,  1970) 

indicated  that  depths  of  the  MED  and  Moho  at  latitude 

5 1  ° 1 1 '  N  and  longitude  96°23'W  are  22.0  and  32.0  km 

respectively.  These  values  are  the  assumed  depths  of  these 
interfaces  at  840  km.  There  is  no  seismic  indication  of  the 
presence  of  the  ARD  in  this  area.  To  conform  with  the 

requirement  of  model  I  the  ARD  is  assumed  to  be  at  a  depth 
of  25.5  km  which  is  the  same  as  that  used  in  the  Swift 
Current-Winnipeg  profile  .  The  densities  used  for  this 

profile  are  the  same  as  the  Swift  Current  -  Winnipeg 

profile,  and  they  are  shown  in  Table  IV.  Other  inversion 
parameters  are  shown  in  Table  4-3. 

Elbow  -  Lake  Winnipeg 

The  original  observed  gravitational  anomalies  is  shown 
in  Fig.  4.23,  and  the  reduced  gravitational  profile  is  shown 
in  Fig.  4.24.  Sampling  interval  was  7  km  and  the  data  series 

contained  258  points. 

Model  Ix  11^  III  and  IV  (Fig.  4.25  to  Fig. 


4.32) 
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The  inverted  shape  of  the  Moho  of  this  profile  differs 
significantly  from  that  of  the  southern  profile.  Besides,  it 
displays  more  relief.  The  crustal  thinning  trend  which  was 
exhibited  by  their  southern  counterparts  is  also  displayed 
by  all  these  4  models.  However,  the  mean  rate  of  crustal- 
thinning  is  much  lower  for  model  III  and  IV;  a  phenomenon 
which  is  also  observed  in  the  Swift  Current  -  Winnipeg 
profile.  The  mean  rate  of  crustal-thinning  for  Model  I  is 
0.014  while  that  of  Model  II  is  0.016. 

The  inverted  depth  of  the  MRD ,  ARD  and  Moho  at  107°50'W 
latitude  are  9.3,  32  and  42.2  km  respectively;  compared  with 
10.4  km  for  the  MRD,  33.4  km  for  ARD  and  43.2  km  for  the 
Moho  at  Swift  Current.  For  Model  II,  the  inverted  depth  of 
the  Moho  is  44.8  km  while  that  of  the  MRD  is  11.9  km.  Like 
their  counterparts  in  the  Swift  Current  -  Winnipeg  profile 
the  inverted  depths  of  the  ARD  of  Model  III,  and  the  Moho  of 
Model  III  and  IV  at  107°50,W  latitude  differ  significantly 
from  those  determined  at  Swift  Current. 
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Table  V 
MODEL  I 


i  - — j - - - - - 1 

|  INVERSION  |  Z  =32. 3km  FRACTIONS.  3  | 

|  PARAMETERS  J  k  =0.087  OFFSET=-250 . 0  mgal  | 

I  I  k  =0.057  ITER= 1 4  | 

|  -H - - - - — - 1 

|  LOCATION  j  ELBOW  j  LAKE  WINNIPEG  | 

|  T  - ( - - - - - - * - - - - ^ 

|  Average  |  A  |  8.38  f  22.0  | 

|  Thickness  of  I  B  |  23.00  |  3.5  | 

|  Layer  (km)  |  C  |  9.80  |  6.5  | 
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MODEL  III 
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INVERSION 
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MODEL  IV 
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INVERSION 

PARAMETERS 


23 .3km  FRACTIONS. 0  | 
=0.087  OFFSET=-240. 0  mgal| 
=0.057  ITER=5  | 


LOCATION 


~r 

I 

I 

-X~ 


I 

■+ 


ELBOW 


|  LAKE  WINN IPEG | 

-1 - i 


Average  Thickness 
of  Layer  (km) 
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Table  V.  Inversion  Parameters  of  Elbow  -  Lake  Winnipeg 
profile. 
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Figure  4.23:  Gravitational  observation  of  Elbow-Lake 
Winnipeg  profile. 


North  of  Winnipeg 


GRAVITY  IN  MGAL 


80 


OBSERVED  GRAVITY 
WTIH  SEDIMENT  REMOVED 


tP 


Figure  4.24:  Reduced  gravitational  observation 
of  the  Elbow-Lake  Winnipeg  profile. 
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Figure  4.25:  Inversion  model  of  Elbow-Lake  Winnipeg 
profile  assuming  model  I.  Vertical  exaggeration  =18. 
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Figure  4.26:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  the  inversion  model  of  Fig  4.25. 
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Figure  4.27:  Inversion  model  of  Elbow-Lake  Winnipeg 
profile  assuming  crustal  model  II. 
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Figure  4.28:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.27. 
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Figure  4.29:  Inversion  model  of  Elbow-Lake  Winnipeg 
profile  assuming  crustal  model  III. 
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Figure  4.30:  Comparsion  between  observed  gravity 
and  gravity  calculated  from  inversion  model  of  Fig  4.31. 
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Figure  4.31:  Inversion  model  of  Elbow-Lake  Winnipeg 
profile  assuming  crustal  model  IV. 
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Figure  4.32:  Comparsion  between  observed  gravity  and 
gravity  calculated  from  inversion  model  of  Fig  4.31. 
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The  southern  extension  of  the  "Nelson  River  Gravity 
High"  corresponds  with  the  ridge  at  445  km 

The  gravitational  anomalies  computed  from  the  4 
inversion  models  agree  well  with  the  observed  gravitational 
observations  iFig.  4.26,  4.28,  4.30  and  4.32). 


Because  of  the 

fact 

that 

at  1 07°50  *  W 

latitude  the 

inverted  depths  of 

the 

Moho 

of  Model  III 

and  IV 

differs 

significantly  from 

those 

determined  at 

Swift 

Current 

(Cumming,  Kanasewich  1966)  Model  III  and  IV  are  not 

considered  acceptable. 

4.5.  Inversion  of  Residual  Gravity  for  Density 

The  preceding  geological  examples  indicated  that,  on 
the  whole,  the  anomalies  computed  from  the  various  models 
fit  the  regional  gravitational  observations  well.  It  was 
also  shown  that  short  wavelength  anomalies  might  probably 
arise  from  the  effect  of  additional  shallow  geological 
features.  If  short  wavelength  components  were  solely 
generated  by  the  shallow  geological  features  how  can  we 
possibly  gain  some  insight  about  the  features  causing  these 
high  wavenumber  components?  As  gravity  interpretation  is 
often  nonunigue  without  any  geological  and  geophysical 
contraints,  there  are  many  possible  features  that  might 
explain  the  same  gravity  anomaly.  This  results  from  the  fact 
that  a  given  gravity  distribution  can  be  produced  by  a 
variety  of  mass  distributions.  For  example,  the  gravity 
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anomaly  that  may  be  explained  by  an  intrusive  body  of 
constant  density  but  variable  thickness  can  also  be  caused 
by  a  similar  body  of  constant  thickness  but  variable 
density.  Thus,  until  the  time  when  additional  information  is 
available,  any  explanation  based  on  the  gravity  data  alone 
is  purely  educated  speculation. 

In  the  following  we  investigate  one  of  the  several 
possible  shallow  features  causing  the  short  wavelength 
components.  A  two-dimensional  slab  of  constant  thickness  is 
presumed  to  be  the  source  of  the  high  wavenumber  components, 
and  we  assumed  variation  of  density  within  the  slab  solely 
responsible  for  the  distribution  of  the  short  wavelength 
components. 

The  Inversion  Scheme 

To  invert  the  difference  between  the  observed  gravity 
and  the  gravity  computed  from  the  inversion  model  for 
density  we  revert  to  equation  (11): 

GO 

F[d  ix)h  (X)  ]  =  -(F[  g  (res)  J/2TTG)  exp  ([  K  ]Z)  -  £  ([  K  ) 

n=2 

F[  d  (X)  hn  (X)  J  (11) 

Where  g  (res)  (residual  gravity)  is  the  difference  between 
the  observed  gravity  and  the  gravity  computed  from  the 
inversion  model.  The  procedure  is  very  similar  to  that 
described  on  page  11  of  Chapter  2,  except  in  this  case  h (x) 
is  presumed  to  be  known  and  the  iterative  scheme  is  used  to 
give  the  density  at  various  locations  along  the  profile. 
There  are  other  minor  differences  as  well.  Since  the 
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densities  obtained  by  inversion  are  of  small  magnitude, 
usually  less  than  2  gm/c.c.  the  convergence  criterion 
(Sn/S1)  for  the  inversion  scheme  has  to  be  considerably 
smaller  than  that  used  to  calculate  h(x).  In  the  following 
inversion  Sn/S2  is  set  <5x10  5.  The  filter  parameters  are 
also  adjusted  to  accommcdate  the  short  wavelength  components 
which  composed  the  major  part  of  the  profile.  Convergence  of 
this  inversion  scheme  is  not  much  affected  by  using  larger 
k  (WH)  or  k(SH)  because  the  residual  gravity  is  usually  of 
small  magnitude;  besides,  there  is  no  power  term  for  d(x). 
Furthermore,  no  offset  is  needed. 

In  the  following  inversions  the  slab  used  for 
investigation  is  5  km  thick  and  the  top  of  this  slab  is  3  km 
below  the  surface. 

Profile  111  1  Swift  Current  -  Winnipeg 

The  residual  gravity,  g (inverted)  -  g  (observed)  of 
Model  I,  and  the  density  obtained  from  inversion  is  shown  in 
Fig.  4.33.  It  can  be  seen  from  Fig.  4.33  that  the  residual 
gravity  and  the  gravity  calculated  from  the  inverted  density 
agrees  with  each  other  very  well.  The  maximum  density 
contrast  is  0.12  gm/c.c.  (at  596  km  of  Figure  4.33).  Thus  at 
596  km,  where  the  local  density  is  2.73  gm/c.c.,  in  order  to 
have  the  gravity  calculated  from  the  inversion  model  to 
agree  with  that  of  the  observed  gravity  the  local  density 
has  to  be  2.85  gm/cc.  It  is  interesting  to  note  that  most  of 
the  activity  is  centered  around  the  southern  extension  of 
the  "Nelson  River  Gravity  High"  (540  km  in  Fig. 4. 33). 
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Figure  4.33:  Density  profile  from  inversion  of  residual 
gravity  of  the  Swift  Current-Winnipeg  profile, and  comparsion 
of  original  residual  aravitv  ^nd  rrravit”  calculated  from 
the  inverted  density  profile. 
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Inversion  parameters  are  shown  in  Figure  4.33. 

Profile  (2) :  Elbow  -  Lake  Winnipeg  (Fig  4.34) 

The  residual  gravity  of  Model  I  of  this  profile  is  used 
for  inversion.  Like  the  preceding  inversion  the  gravity 
calculated  from  the  inverted  densities  fits  the  residual 
gravity  very  well.  Two  major  density  anomalies  occur  in  this 
profile.  The  first  one  centres  around  485  km,  and  is  in  the 
region  of  the  southern  extension  of  the  "Nelson  River 
Gravity  High".  The  second  one  is  around  800  km,  just  west  of 
Lake  Winnipeg.  The  significance  of  these  density  anomalies 
is  discussed  in  the  following  section. 

4.6.  Interpretation 

As  seen  from  the  cross  sections  of  the  inversion  models 
of  profile  (1)  and  (2),  the  inverted  topography  of  the  Moho 
of  Model  I  shows  best  agreement  with  that  obtained  from 
seismic  data.  The  inverted  topography  of  the  Moho  of  Model 
II  is  very  similar  to  that  of  Model  I,  however,  the  mean 
slope  of  the  Moho  tends  to  be  greater  for  Model  II. 
Structural  relief  is  considerably  stronger  in  Model  III  and 
IV;  a  phenomenon  which  is  consistent  in  both  profile  (1)  and 
(2)  . 


Gravity  analysis  of  the  Swift  Current  -  Winnipeg 
profile  reveals  theat  the  thickness  of  the  crust  decreases 
eastward  at  a  mean  rate  of  0.013.  Other  than  the  depression 
at  the  west  end  of  the  profile,  the  relief  of  the  Moho  in 
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Figure  4.34:  Density  profile  from  inversion  of  residual 
gravity  of  the.  Elbow-Lake  Winnipeg  profile, and  comparsion  of 
original  residual  gravity  and  gravity  calculated  from  the 
inverted  density  profile. 
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the  interior  of  the  Churchill  province  is  flat.  The  Moho 
exhibits  little  relief  in  the  Superior  province,  except  for 
the  depression  at  the  east  end  of  the  profile.  The  southern 
extension  of  the  "Nelson  River  Bouguer  gravity  High”  is 
inverted  as  a  structural  high.  The  depressions  of  either 
side  of  this  structural  high  may  be  due  to  down  buckling  of 
the  crust.  The  same  might  have  happened  for  the  depressions 
at  either  end  of  the  profile.  Bell  (1971b)  believed  that  the 
surfacing  of  the  Pikwitonei  province  was  probably  due  to  the 
up-tilting  of  the  continental  plate  which  was  made  up  of  the 
Pikwitonei  basement  complex  and  the  Archean  supracrustal 
rocks  of  the  Superior  province.  The  inverted  structural  high 
seems  to  be  in  accordance  with  Bell's  (1971b)  idea. 

The  inverted  topography  of  the  Moho  of  the  Elbow  -  Lake 
Winnipeg  profile  differs  significantly  from  that  of  the 
Swift  Current  -  Winnipeg  profile.  The  structural  high  which 
occurs  in  profile  (1)  is  also  present  in  this  profile;  but 
of  lesser  relief.  West  of  this  topographic  high,  the  Moho 
displays  much  more  relief  than  that  of  profile  (1).  However, 
east  of  this  structural  high  the  relief  of  the  Moho 
resembles  its  southern  counterpart.  The  mean  rate  of  crustal 
thinning  is  similar  to  that  of  the  Swift  Current  -  Winnipeg 
profile. 

Looking  at  the  inverted  density  waveform  of  the  Swift 
Current  -  Winnipeg  profile  and  the  Elbow  -  Lake  Winnipeg 
profile  one  can  observed  that  all  the  high  amplitude 
components  occur  in  the  eastern  portion  of  both  profiles. 


, 

■ 
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This  is  especially  obvious  on  profile  (2).  (Fig.  4.34)  In 
fact  both  profiles  can  be  roughly  divided  into  two  sections: 
the  eastern  section  with  the  high  amplitude  density 
contrast,  and  the  western  section  with  the  low  amplitude 
density  contrast.  The  dividing  line  for  profile  (1)  (Fig. 
4.33)  is  at  432  km,  while  that  of  profile  (2)  was  chosen  at 
370  km  .  The  range  of  density  contrast  of  profile  (1)  is  - 
0.09  gm/cc  to  +0.12  gm/cc,  which  gives  a  density  range  of 
2.64  to  2.85  gm/cc  .  On  the  other  hand,  in  profile  (2)  the 
range  of  density  contrast  is  -0.13  to  0.23  gm/cc,  and  its 
density  ranges  from  2.6  gm/cc  to  2.96  gm/cc  .  What  possible 
information  can  be  gathered  from  these  observations?  Do  all 
these  density  anomalies  arise  from  lithological  contrast 
within  the  basement,  or  are  they  pseudo-density  anomalies 
reflecting  structural  variation  within  the  basement?  Little 
information  concerning  basement  geology  and  structure  of  the 
area  is  available,  and  it  is  not  possible  to  state  whether 
the  gravity  anomalies  giving  rise  to  the  negative  density 
contrasts  are  mainly  due  to  low  density  rocks  like  schist 
(2.30-2.87  gm/c.c.),  granulite  (2.57-2.73  gm/c.c.),  gneiss 
(2.59-3.0  gm/c.c.),  or  are  mainly  structural  features  (like 
faults,  shear  zones  and  grabens)  •  If  negative  density 
contrasts  were  caused  by  low  density  rocks  these  density 
anomalies  will  have  lithological  significance,  however,  if 
they  originated  from  faults  and  such  like  they  will  have 
structural  significance.  Likewise,  it  is  difficult  to  say 
whether  the  positive  density  contrasts  originated  from  dense 
rocks  (gabbro,  diorite,  dykes,  sills  etc)  or  from  structural 
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features  (steps,  horsts  ect) .  It  is  most  likely  that 
lithological  and  structural  contrasts  within  the  basement 
are  both  responsible.  In  his  delineation  of  the  boundary 
between  the  Pikwitonei  province  and  the  Superior  province 
Bell  (1971b)  used  the  distribution  of  mafic  dyke  swarms  to 
substantiate  his  interpretation  of  the  boundary.  Several 
distinct  dyke  complexes  were  identified  or  inferred.  The 
north-northeast  gabbro- diabase  dykes,  called  the  Molson  Dyke 
Swarm  (Fabrig,  et  al. ,  1965,  P289) ,  occur  in  large  numbers 
in  the  rocks  at  the  west  end  of  Cross  Lake  subprovince  (S.W. 
of  Pikwitonei  province;  belongs  to  the  Superior  province) 
and  the  granulites  of  the  Pikwitonei  province.  However,  none 
of  these  dykes  are  known  to  intrude  any  of  the  Churchill 
province2  rock  lying  northwest  and  north  of  the  Pikwitonei 
complex  gneisses  (Bell  1971b) .  This  fact  was  utilized  by 
Bell  (1971b)  to  extend  the  Superior-Churchill  boundary 
beneath  the  Paleozoic  cover  lying  west  of  Lake  Winnipeg. 
Bell  (1971b)  also  pointed  out  the  almost  exact  parallelism 
between  the  Churchill-Superior  boundary  and  the  'Molson 
Dykes  ' . 

Two  distinct  dyke  complexes  are  also  aligned  with  the 
Pikwitonei-Churchill  boundary:  Cuthbert  and  Nelson  dykes 
(McDonald,  1960;  Patterson,  1963,  pp  27-29).  They  grade  from 
peridotite  to  mafic  gabbro,  and  intrude  rocks  of  the 
Pikwitonei  province  and  the  Cross  Lake  subprovince  (Bell 


2  Bell's 
subprovince. 


Churchill 


province 


includes  the 
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1971b)  . 

Bell  (1971b)  postulated  that  ultramafic  dykes  (?)  are 
responsible  for  the  coincidence  of  the  southwestern 
extension  of  the  'Nelson  River  Bouguer  Gravity  High'  with 
the  aeromagnetic  highs  at  longitude  100°00'W  and  latitude 
53°30'N.  Bell  (1971b)  called  it  the  Cedar  Lake  dyke(?).  Bell 
(1971b)  also  postulated  that  the  Pelican  Lake  dyke  (?)  , 
which  is  50  km  to  the  southeast  of  the  Cedar  Lake  dyke,  is 
responsible  for  the  chain  of  linear  magnetic  anomalies, 
which  is  fringed  with  a  gravity  high  to  the  northwest.  Bell 
(1971b)  suggested  the  Pelican  Lake  dyke  (?)  is  similar  to  the 
Fox  River  ultramafic  sill  or  sills,  and  the  Pelican  Lake 
dyke  is  just  west  of  the  Churchill-Superior  boundary. 

As  already  indicated  there  are  a  number  of  high- 
amplitude-positive-density  contrasts  around  the  southern 
extension  of  the  'Nelson  River  Bouguer  Gravity  High'.  These 
density  anomalies  probably  correspond  with  dense  rocks, 
dykes  or  sills.  The  density  anomalies  east  of  600  km  of  both 
profiles  most  likely  belong  to  the  Superior  province, 
because  they  are  in  an  area  with  E-W  trends;  a  trend  typical 
of  the  Superior  province.  In  area  A  of  profile  (1)  and  area 
(B)  of  profile  (2)  there  are  observable  geophysical 
resemblances  to  that  displayed  by  the  Cedar  Lake  dyke(?)  and 
the  Pelican  Lake  dyke.  Positive  density  anomalies  of  both 
areas  A  and  B  correspond  with  gravity  highs  which  either 
have  coincident  magnetic  highs  or  fringing  magnetic  highs. 
(Federal-Provincial  aeromagnetic  map  coverage  ends  near 
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latitude  52°N  and  observation  was  made  partly  from  the 
magnetic  map  in  appendix  C  and  partly  from  a  confidential 
map  of  the  Gulf  Research  and  Development  Company) .  If  the 
western  limit  of  the  Superior  province  at  latitude  52°30*N 
and  longitude  100°W  is  terminated  byN  the  Pelican  Lake 
dyke  (?)  (Bell  1971b)  point  C  and  D  of  profile  (1)  and  (2) 
respectively  probably  mark  the  western  limit  of  the  Superior 
province  as  well. 

Like  the  magnetic  and  gravity  anomalies  of  the 
Churchill  province,  which  are  small  and  low  amplitude  the 
density  profile  west  of  point  E  of  profile  (1)  and  F  of 
profile  (2)  does  not  show  the  high  amplitude  density 
anomalies  which  are  displayed  in  the  eastern  portions.  It 
was  indicated  that  dykes  and  sills  are  abundant  in  the 
Pikwitonei  province,  and  also  that  some  of  the 
metasediments,  meta volcanics ,  and  migmatic  gneisses  of  the 
Waboden  subprovince  are  of  lower  to  upper  amphibolite 
facies.  The  high-amplitude-positive-density  anomalies  in 
area  A  of  profile  (1)  and  B  of  profile  (2)  afford  the 
possibility  of  similar  occurrence  of  these  dense  rocks. 

The  dividing  line  between  the  Churchill  province  and 
the  boundary  zone  is  postulated  to  be  at  point  E  of  profile 
(1)  and  F  for  profile  (2) .  For  Superior  province  we 
postulate  points  C  and  D  of  profile  (1)  and  (2)  respectively 
to  be  the  western  limit  of  its  existence. 
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4.7,  Conclusions 

The  inversion  scheme  developed  by  Parker  and  Oldenburg 
(1974)  has  been  extended  to  cover  a  multi-layered  medium. 
Unlike  other  conventional  gravity  inversion  methods  this 
multilayer  inversion  scheme  not  only  takes  into  account  the 
gravitational  effects  of  strata  overlying  the  inversion 
interface  but  also  allows  for  lateral  variation  of  density 
and  thickness  within  each  layer.  Interfaces  overlying  the 
inversion  horizon  are  assumed  to  have  the  same  shape  as  the 
inversion  interface.  It  has  been  shown  that  if  variation  of 
density  and  thickness  of  beds  overlying  the  inversion 
interface  is  only  known  at  certain  locations  of  the  profile 
a  linear  interpolation  of  these  parameters  still  gives 
satisfactory  inversion  models.  Inversions  based  on  various 
crustal  models  indicate  that  the  inverted  topography  of  the 
inversion  interface  of  inversions  models  which  take  into 
account  of  the  effect  of  overlying  layers  often  shows  lesser 
relief  than  that  obtained  from  inversion  models  which 
attribute  all  gravitational  variation  to  a  single  surface. 

The  nonuniqueness  of  the  gravity  inversion  is 
characterised  by  two  free  parameters:  d,  the  density 
contrast  between  adjacent  layers,  and  Z  ,  the  level  at  which 
the  inversion  is  made.  Without  additional  information 
constraining  these  two  parameters,  the  ambiguity  in  the 
gravity  interpretation  cannot  be  reduced. 


The  Parker-Oldenburg  (1972,  1974)  algorithm  has  been 
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rearranged  to  find  the  density  distribution  of  a  near 
surface  slab  in  order  to  fit  the  difference  between  the 
original  gravity  and  the  gravity  computed  from  the  inverted 
topography.  Together  with  other  geological/geophysical 
controls  the  resultant  density  waveform  provides  a  valuable 
tool  for  interpretation.  However,  in  the  absence  of  other 
information  caution  must  be  exercise  in  interpretating  the 
resultant  density  profile.  This  is  because  a  given  gravity 
distribution  can  be  produced  by  a  variety  of  mass 
distributions. 

The  four-layered  inversion  model  of  the  Swift  Current  - 
Suffield  profile  is  in  reasonable  agreement  with  the  crustal 
model  proposed  by  Cumming  and  Kanasewich  et  al  (1966). 
Gravity  analysis  of  the  Swift  Current  -  Winnipeg  profile  and 
the  Elbow  -  Lake  Winnipeg  profile  reveals  a  crustal  thinning 
towards  the  east.  The  southern  extension  of  the  "Nelson 
River  Bouguer  Gravity  High”  is  inverted  as  a  structural  high 
which  probably  corresponds  with  the  up-tilting  of  the 
Pikwitonei  province  (Bell  1971b). 

Density  analysis  of  the  residual  gravity  of  the  Swift 
Current  -  Winnipeg  profile  and  the  Elbow  -  Lake  Winnipeg 
profile  reveals  several  high  amplitude  anomalies  in  the 
central  and  eastern  portions  of  both  profiles.  Some  of  the 
positive  density  anomalies  around  the  location  of  the 
southern  extension  of  the  "Nelson  River  Bouguer  Gravity 
High”  probably  correspond  with  dykes/sills  which  are  so 
abundant  in  the  Pikwitonei  province,  or  with  the  amphibolite 
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facies  of  the  Wabowden  subprovince.  For  the  Swift  Current  - 

\ 

Winnipeg  profile  the  western  and  eastern  limits  of  the 
boundary  zone  are  tentatively  set  at  102°15*W  longitude 
50°20'N  and  100°00  longitude,  50°20*  N  respectively.  As  for 
the  Elbow  -  Lake  Winnipeg  profile  the  eastern  limit  of  the 
Churchill  province  is  set  at  longitude  102°35,W  lat. 
51°10.5'N  and  the  western  limit  of  the  Superior  province  is 
set  at  longitude  101°W  latitude  51°10.5'N  (Figure  4.35). 


Figure  4.35:  Gravity  map  of  the  southern  plains  of 
Saskatchewan  and  Manitoba  showing  possible  locations  of 
the  western  limit  of  the  Superior  province  and  the  eastern 
limit  of  the  Churchill  province. 
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APPENDIX  A 

Derivation  of  Parker-Oldenberg  Algorithms 

In  the  following  discussion  we  denote  any  three- 
dimensional  position  vector  with  a  bar  (e.g.  "r  )  and  its 
projection  on  the  X-Y  plane  by  an  arrow  (e.g.  ~r  )  . 

The  upper  and  lower  boundaries  of  a  single  layer  of 
material,  which  is  contained  in  a  finite  area  D,  is 
described  by  Z=h(~r)  and  Z=0  (a  flat  surface)  respectively. 

He  further  assume  h  fr)  is  bound  and  integrable.  We  also 
define  the  two-dimensional  Fourier  transform  of  any  function 
g  (r)  and  its  inverse  by: 

F[  g  (T)  J  =f  dA  g  (rj  exp  (iK^T) 

*x 

g(T)=(1/47 TZ)J^  dw  F[g(r)  ]  exp  (-iT-T) 

respectively.  X  indicates  integration  over  the  X-Y  plane  and 
K  integration  over  the  Kx~Ky  plane,  dA  =  dxdy,  and 
dW=dKx  dKy. 

Consider  a  point  on  a  horizontal  plane  which  is 
everywhere  above  h(r)  and  is  at  a  distance  Zo  above  the 
origin,  the  potential  field  at  a  position  r0  from  the  origin 
is  given  by: 
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U(rQ)  =  Gd  /  dV/  |  rG  -  r  j 


=  Gd 


/»  /»z=h(r) 

/  dA  j  dz/| rn -r | ; 

«/D  Jz-0 


(A.  1) 


where  G  is  the  Gravitation  constant;  d  is  the  density  of  the 
layer,  which  is  assumed  constant;  7  is  the  position  vector 
from  the  origin  to  the  source;  and  D  is  the  finite  area  in 
the  X-Y  plane. 


The  two-dimensional  Fourier  transform  of  equation  (A.1) 


is 


F[U(r0)j  =[  dA0U  (f0  )  exp  iiT*  r  0  ) 

x 


■/, 


=  Gd  J  dA0  exp liK- r0) 


f  dAf  ” 

Jo  J2=0 


h(?T 


dZ/| rft -r | 


Note:  TT«r  =  K«rt 


Interchange  the  order  of  integration  we  obtain: 


F[U]=Gd  fQ 


exp  (iK «  ro )  ]/|ro-r|  . 


(A. 2) 


The  last  integral  of  equation  (A. 2)  can  be  expressed  as: 


exp(iT-r£)  ]/|r0-r|  =  exp(iT-T) 
exp(i|T|  |To-T|cos@  )  J/l  r0 -  r  | 


A  A 

r0  =  r0-Z(Z-r0) 

=  to  ~  ( Z  •  zj  ro  « 


(A. 3) 


because : 
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Q  is  the  angle  between  the  two  vectors  IT’  and  fro - r"j  in  the 
X-Y  plane.  Let  u  =  I'rJj-Irl  and  transform  (A. 3)  to  polar  co¬ 
ordinates  we  have: 


_ r00  _ 

=  exp(iK*r)/  d$/[udu  exp(i|K|ucos  e 

Jo  Jo  Q 


—  —  J2W  SO 

expiiK*r )j  dQj  [  udu  exp  (i  j  K  l  ucos  Q  )  ]/[  u2  +  (ZQ-Z)  J 


i/2 


If  we  let  B=Z0-Z,  the  integral  becomes: 


— *._»  . /p /*2"^  — » 

=  expiiK*r)  /  udu/iu2+B  2)'l  [  exp  (i  1  K  [  ucos  Q  Jd  Q  (A.  4) 

** 0  •'o 


Since  J0  ( I K | U  ),  which  is  a  Bessel  Function,  is: 


=  1/2  77' [  exp  (i  |T|  ucos  8  ]  d  Q  , 


(A. 4)  now  becomes: 


=  2  7 r  exp  (iK 


00 

’~r)  /  udu 
Jn 


J0(|K|  u)  /(u2  +  B2) 


1/2 


(A. 5) 


The  integral  of  (A. 5)  can  be  found  from  standard  integral 
tables,  and  is: 


-  [  2 -nr  /|F|  ][  exp(iK-T-|K|  (Z0-Z)  J 


Substituting  (A. 6)  back  into  (A. 2)  we  have: 


r  fir) 

ll  dA  /  1 

Jo  Jo 


F[  U  J  =  2TTG  d dA  J  dZ[  exp  (iK*  r- I  K  |  (Z0-Z)  J/|  K  | 


(A.b) 


l  A  .  7 ) 


. 


no 


Integrating  over  Z  from  0  to  h  (?)  equation  (A. 7)  becomes: 


•L 


F[U]  =  27TGdJD  dA[  exp(iK.r-jK|Z0)  ][exp(Kh(r)-1)  ]/|K|2 


(A. 8) 


—  —  00  __ 

Since  expiKh(r)-l)  can  be  expressed  as  £  [  I  K | n- 2hn  (r)  ]/n ! 

n-  i 

in  Taylor  series  equation  (A. 8)  can  now  be  expressed  as: 


F[UJ  = 


—  /»  _^_*00  n 
27TG  d[expi-|K|Z  ]|  dAdexp  (iK-  r)  V  [  |K|n-2hn  (?)  J/n! 

Jo  4 


Since  | K |  is  a  constant  with  respect  to  dA: 


F[U]  =  277Gd  [  exp-  jlT  |  Z0  ]X  |”1T  |  n  —  2 /n  !  /*  dAhn  (?)  exp  (iK*  r) 

n=|  o'D 


Integrating,  we  obtain  equation  (A. 9)  which  is: 


F[  0  ]  =  2 7713 d  exp  (-|T|  Zft)  Y  |~K  |  n- 2/n !  F[  hn  (r)  ] 

0  n=l 


(A. 9) 


To  find  the  vertical  attraction  of  the  perturbing  mass, 
we  note  that  for  any  observation  datum  above  the 
masses  y 2U  =  0.  This  allows  us  to  write: 


U 


<r0)  =  l/mrzj*  fl2KIJ  (?)  exp  (-iK-^  ) 


To  upward  continue  the  above  expression  to  the  observation 

— «  A  _ 

datum  we  have  to  multiply  it  by  exp  (- | K | Z* r0) . 


0  (r0 )  =  1/4TT2  /  d2 KU  (?)  exp  i-iK  *ro-i|  K|  Z*  r0  ) 


. 


Ill 


Thus  F[U(^)J  =  U  (T)  exp  i-|T|  Z- r0)  . 

By  definition,  the  vertical  attraction,  Ag,  is: 

Ag=+  du/  d  z. 

=  F[AgJ  =  -mnuj 

_  00  _ 

=  F[i  g]  =  — 27TGd  [exp  1-1  K|150)  ]  X  |  K  |  n“  »  F[  hn  [r)  J/n ! 

n=l 

which  is  Parker* s  algorithm. 


■ 
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APPENDIX  B 


Convergence  of  Parke r-Qldenburg  Algorithms 
Parker  (1972)  Algorithm: 

We  revert  to  Parker's  (1972)  algorithm: 

—  —  00  n  _ 

F[  A  g  (r)  J=~27TGd  exp  (-  |  K|  Zo)  £  |  K  j  n- i/n !  F[  hn  (r )  J  ;  (B.1) 

n=| 

Mr)  is  assumed  to  be  bounded  and  integrable,  and  vanishes 
outside  some  finite  domain  D.  Two-dimensional  Fourier 
transform  of  h(r)  is: 


F[  hn  (r)  ]  =  f 
J  D 


dAhn  fr)  exp  (iK  •If)  ,  and: 


dA | hn (r)  | «  j  exp (iK*T)  | 


d  A  J  hn  i~r)  | 


<  AH 


n 


A  =  area  of  D 

H  =  maximum  value  of  Jh(r)  | 
Substituting  the  bounds  into  equation  (B.1)  we  have: 

|F[Agir)]l  S  |-2W5d  /fKlJ  (I  K  |  H)n  A/n !  | 

n«l 


'V  * 
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I  Ft  A  g  IT)  ]|  <  27TGd  A/m[exp(|F|H)-1  J;  (B.2) 

Equation  (B.2)  is  uniformly  and  absolutely  convergent 
in  any  bounded  domain  of  the  K-plane  (Whittaker  and  Watson, 
1962,  p581 ) . 

To  show  the  rate  of  convergence  can  be  optimised  we  let 

00 

S  =  £  (K  n/n  I )  exp  (-  KZo)  F[  hn  fr)  J  and  write  K  for  |K|.  Since 

nsi* 

| F[  hn  (r)  ]|  1  AHn  we  have: 

00  «  - 

S'  =  Aexp  (-KZo)£  k  H/n!' 

n=i 

and  | S |  <  S'. 

The  preceeding  result  can  be  shown  easily.  If  we  take 
the  nth  term  of  both  functions  we  have: 

Sn  Kn/n  !  F[  hn  (T)  J 

S'n  A  Kn/n!  Hn 

Sn  F[  hn  (T)  ] 

S'  n  AHn 

<1  V  |  F[  hn  fr)  J  |  <  AH° 

There  is  an  upperbound  for  S*  .  The  nth  term  of  S'  is: 

Sp  =  [  A  (KH)n /n!  ]/exp  (KZo) 

00 

=  [  A  (KH)n/n!  ]/  £  (KZo)'  /i ! 

i=0 

<  [  A  (KH)n  /nl  ]/[  (KZo)n  /n!  ] 

<  AHn/zJ 

<  A(H/Zo)R 


■ 
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Summing  over  all  terms  we  have: 

00 

S'  <7  A  (H/Zo)n 
n=l 

Therefore  when  H<Zo  and  Zo>0  the  serie  for  S'  is 
uniformly  convergent  in  the  whole  K-plane,  as  demonstrated 
by  Weierstrass  M-test  (Whittaker  and  Watson,  1962,  p49)  and 
hence  this  is  true  for  S.  The  smaller  H/Zo  can  be  made,  the 
faster  will  be  the  rate  of  convergence. 

It  is  not  apparent  that  we  have  any  control  over  H/Zo 
in  a  given  calculation,  but  this  is  in  fact  not  the  case.  In 
choosing  our  coordinate  system,  we  chose  Z=Q  to  be  the 
bottom  of  the  layer  of  material.  It  is  an  arbitrary  choice, 
for  we  may  choose  any  depth  for  an  origin.  Deplacing  the 
origin  does  not  affect  the  validity  of  equation  (B. 1)  but  it 
does  change  the  numerical  values  of  Zo  and  h  (r) ,  and  also  H. 

Oldenburg  (1974)  Algorithm: 

Similar  approach  was  used  to  derived  at  the  convergence 
criterion  of  Oldenburg's  inversion  scheme. 

We  revert  to  Oldenburg's  inverse  formula: 

F[  h  (X)  ]  =  -F[Ag(X)  Jexpd  K|Zo)/27TGd 

-V  (I  K|  n-*/n! )  F[  hn  (X)  ];  (B  .3) 

n=2 

The  one-dimensional  Fourier  transform  of  a  function  h(x) 


by: 


' 


■ 
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F[hn(x)]  =[  dxhn  (x)  exp  (iKx) 

J  D 


and  |  F[  hn  (x)  ]  j  =  if  dxhn(x)exp(iKx)| 


Jo 

<  /  lhn 

•'D 


(x)  |  dx 


<  L  H 


n 


L  is  the  length  of  D  and  H  =  max|h(x)  |.  We  also  require  that 
hix)  be  bounded  and  integrable  and  h(x)=Ofor  X  and  D. 

Inserting  these  bounds  into  the  infinite  sum  in 
equation  (B.3)  we  get: 

IV  iKH-i/n!)  F[hn  (X)  J|  <  (L/K)  £  KnHn/nl 
2  n=2 

=  iL/Kj  [  exp(KH) -1-KH J;  lB.4) 


The  right  hand  side  of  iB.4)  is  bounded  for  any  finite 
value  of  K,  hence,  from  the  properties  of  the  exponential 
function  (Whittaker  and  Watson,  1962,  p581)  the  sum  of 
Fourier  transforms  is  absolutely  and  uniformly  convergent  in 
any  bounded  domain  of  the  K-plane. 

The  convergence  criterion  chosen  by  Oldenburg  is  that 
successive  terms  in  the  sum  are  computed  until  Sn/S2  <  5E-3. 
For  the  forward  algorithm,  sufficient  number  of  terms  are 
evaluated  untill  Rn/RI  <  10E-6. 


. 


? 


'  I 
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max 

Sn  =  over  all  |  ( |  K|  i/n  ! )  F[  hn  (x)  ]  | 
K 


and 


max 


Rn  =  over  all  |  exp  (- |  K  |  Zo)  |  K|  n_1/nl  F[  hn(x)  J| 
K 


Though  convergence  of  the  infinite  sum  is  assured, 
however,  there  is  still  no  control  over  the  convergence  of 
the  iterative  procedure.  Oldenburg  was  unable  to  establish 
analytical  convergence  criterion  for  the  iterative 
procedure.  He  resorted  to  an  empirical  criterion  which 
requires  that  the  mean- root- square  difference  between  two 
consecutive  approximations  of  h (x)  be  less  than  some 
arbitrary  chosen  small  number.  The  mean- root-square  error 
can  be  expressed  mathematically  as: 


Error 
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Magnetic  map  of  southern  Saskatchewan. 


MAGNETIC  ANOMALY  MAP 
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Computer  program  to  do  multilayer  gravity  inversion. 


C  ***#*******♦♦##♦«**#*#♦***#**** **#**#*# 

C  MULTILAYER  GRAVITY  INVERSION  USING  FOURIER  TRANSFORMS 

C  PROGRAM  MODIFIED  FROM  OLDENBURG'S  INVERSION  PROGRAM 

C  WRITEN  BY  J.S.K.  LEE;  UNIVERSITY  OF  ALBERTA,  1977. 

C  PROGRAM  TO  CALCULATE  THE  GR AVITATINN AL  ANOMALY  AND  DO  THE 

C  INVERSION  OF  MULTILAYER  MODEL,  WITH  LATERAL  LINEAR 

C  VARIATION  OF  DENSITY  AND  BED  THICKNESS. 

C  INPUT  TO  THE  PROGRAM  SHOULD  BE  OF  THE  FOLLOWING  FORMAT 
C  **  N2,NN,IG,NJJ,NLL,NKK,NLAY 

C  (#  OF  OBSERVED  GRAVITY,  #  OF  DIGITISED  STATIONS  TO  BE  PRODUCED, 
C  IG  .NE.  1  FOR  DATA  NOT  IN  TENTH  MGAL ,  #  OF  FILTERS  TO  BE  USED, 
C  #  OF  DENSITY  VALUES  FOR  INVERSION  HORIZON ,  #  OF  Z,  AND  TOTAL 
C  #  OF  LAYERS. 

C  **  BTK  (BEDTHICKNESS  OF  OVERLYING  BEDS  AT  W  END, FROM  THE 
C  INVERSION  HORIZON  UP, TOTAL  INPUT  SHOULD  =  NLAY-2  ) 

C  **  BTK2  (SIMILAR  TO  BTK  BUT  AT  THE  E  END) 

C  **  DERHO  (DENSITY  OF  OVERLYING  LAYERS  AT  W  END, IN  ASCENDING 
C  ORDER  FROM  THE  INVERSION  HORIZON) 

C  **  DERH02  (SIMILAR  TO  DERHO  BUT  AT  E  END) 

C  **  FRAC, UM  (FRACTION  OF  INPUT  GRAVITY  TO  BE  USED  FOR  FIRST  ITER 
C  AND  AMOUNT  OF  OFFSET  TO  BE  USED  TO  FACILIATE  CONVERGENCE) 

C  **  ZUP  (INVERSION  HORIZON,  NKK  TERMS) 

C  **  RHO  (DENSITY  CONTRAST  OF  INVERSION  HORZN,  NLL  TERMS) 

C  **  RHOTWO  (DENSITY  CONTRAST  OF  INVSN  HORZ  AT  E  END,  NLL  TERMS) 

C  **  JPRINT  (INDEX  FOR  PRINTING  INTERMEDIATE  RESULTS,  1=YES,  3=NO 
C  **  X(J)  (X  CO-ORDINATE  OF  INPUT  GRAVITY) 

C  **  GG  (J)  (OBSERVED  GRAVITY) 

C  **  NR  (#  OF  SEDI-ENT  GRAVITY) 

C  **  GG  (J)  (SEDIMENT  GRAVITY  CALCULATED  AT  THE  SAME  LOACATIONS  AS 
C  OBSERVED  GRAVITY) 

C  **  W rl , S H  (LOW  CUT-OFF  AND  HI  CUT-OFF  WAVENUMBERS,  NJJ  TERMS) 

Q  **$*$**************************************♦*♦*** 

DIMENSION  H  1  (5  12)  ,  H  (512)  ,  X  (51  2)  ,  GR  A  V  (5 1  3)  ,DG  (512) 

DIMENSION  DG3  (51 2)  , BTK  (2 0) , DG2  (512)  ,G1  (512)  ,G2  (512) 
DIMENSION  RHO  (20) , BTK2  (2  0) , DERHO (2C) , DG4 (512) ,ZZ (2  0) 
DIMENSION  RODIG1  (512)  ,RODIG2(512)  ,F(512),H2(512),B(512) 
DIMENSION  DERH02 (20) , RHOTWO (20) 

DIMENSION  AA  (512)  , BB  (512) ,CC (512)  ,GG  (512) 

COMPLEX  SUM  (512)  ,  HTON  (512)  ,  ANOM(512) 

COMMON/BOUG/ ACC, IPRINT 
COMMON/HI LLS/NN,DX, ZUP, GRHO,GRO 

COMMON/IN V/ITE RMS, ITER, ISKIP, AC 1 ,IZUP, JPRINT, FRAC, NLA Y 
DATA  PI/3.1415926535/ 

DATA  G/6.67/ 

READ (5, 100)  N2,NN,IG,NJJ,NLL,NKK, NLA Y 
100  FORMAT (714) 

NBD=NL AY-2 


. 


on  no 
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(ZZ  (I)  ,1  =  1  f  NKK) 

(RHO (I) ,1=1 , NLL) 
(RHOTWO  (I)  ,1=1 , NLL) 
JPRINT 


READ  (5,  101)  (BTK(T)  ,1=1, NBD) 

READ(5,101)  (BTK2 (I) ,1=1, NBD) 

R EAD (5 , 101)  (DERHO  (I)  ,1=1, NBD) 

READ  (5, 101)  (D  ER H02  (I)  ,1  =  1,  NBD) 

101  FORMAT  (10F8. 3) 

READ  1 5 , 1 7 1 )  FR  AC , UM 

171  FORMAT (2®. 3) 

READ  (5,102) 

READ (5,102) 

READ (5,102) 

READ  (5,222) 

IPRINT=J PRINT 

222  FORMAT  (II) 

102  FORMAT  (5F8. 3) 

:  INPUT  THE  #  OF  TERMS  TO  BE  USED  IN  POWER  SERIES 

NTERMS-35 

READ  (5,47)  (X  ( J)  ,  J=  1 ,  N2) 

47  FORMAT  (10F8. 3) 

READ  (5,4  7)  (GG  (J)  ,J=1,N2) 

:  IF  INPUT  GRAVITY  IS  TENTH  MILLIGAL  CONVERT  TO  MILLIGAL 

IF  (  G.EQ. 1)  GO  TO  99 
DO  1  1=1, N2 
1  GG (I) =GG  (I) *0.1 

11  FORMAT  ( 10F1  2.  1) 

99  WRITE  (6,56)  ( X ( J)  , GG  ( J)  ,  J=  1 ,  N2) 

56  FORMAT  (/////, 40X , 1  INITIAL  DISTANCE  IN  KM.  AND  GRAVITY  ANOM 
1  //, 10 (2F9.3,6X,  2F9.3,6X ,2F9. 3 , 6X, 2F9 . 3 , 6X, 2F9 . 3,/) ) 

309  FORMAT (/////, 30X, ’  SEDIMENT  GRA VIT Y * , //, 1 0 ( 10 F 1 2. 

1  **,/)) 

:  CALCULATE  DIGITISING  INTERVAL 

XM AX=X (N2) 

DX=X  MA X/NN 

##****#****##***#**#***************#*# 


CALL  SPLNCO  TO  PRODUCESPLINING  COEFFICIENTS. 

CALL  S  PLNCO (1,N2,X, GG, AA,BB,CC) 

DO  23  J= 1 , NN 
H  ( J) = (J-1)  * DX 
11  =  0 

CALL  SPLNEV  TO  CALCULATE  SPLINED  VALUES  AT  EACH  POINT. 
SPLINED  VALUES  ARE  STORED  IN  VECTOR  G2. 

CALL  SPLNEV  (1,N2, X, GG, AA, BB , CC , H ( J )  ,G2  (J)  ,11) 

23  GRAV  ( J)  =G2  (J) 

WRITE  (6,55)  (H(J)  ,GRAV(J)  ,J=1,NN) 

C  INPUT  SEDIMENT  GRAVITY  AND  DIGITISE 
READ (5, 105)  NR 
105  FORM  AT (14) 

READ  (5,4  7)  (GG  ( J)  ,  J=  1  ,  N R) 

IF  (JPRINT.LT.  3)  WRITE  (6,309)  (GG  ( J)  ,  J=  1  , N R) 

CALL  SPLNCO (1 , NR, X , GG, A A , BB , CC ) 

DO  25  J=  1 , NN 
H  (J)  =  (J-1)  *DX 
11=0 

C  CALL  SPLNEV  TO  CALCULATE  SPLINED  VALUES  AT  EACH  POINT. 


. 
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C  SPLINED  VALDES  ARE  STORED  IN  VECTOR  G1. 

CALL  SPLNEV (1,NR,X,GG,AA,BB,CC,H(J)  ,  G  1  (  J)  ,11) 

25  GRAV(J)  =GRAV  (J) +G1  (J) 

IF  i  JPRINT.  GE.  3)  WRITE  (6,905)  (H  ( J)  ,  J=1  ,  NN) 

905  FORMAT (///// ,  30X, *  DEGITISING  INTERVAL® , //, 10(10F12. 

1  <*,/)) 

IF  (JPRINT.  GE.  3)  WRITE  (6,509)  (GRA  V  (J)  ,  J=  1  ,  NN) 

501  FORMAT  (/////, 40  X, *  DISTANCE  IN  KM.  AND  SYM.  GRAVITY  IN  MILL 
1 r//r 10 (2F9.3,6X,  2F9.3,6X,2F9.3,6X,2F9.3,6X,2F9.3,/) ) 

C  *#*********#**##*##*##****###***#:$:#######  **#*#«*## 

C  TO  GENERATE  A  SYMMETRICAL  SERIES. 

C  ♦♦♦♦A******##*###*****#*******#*######*#*##*#** 

DO  71  J=1,NN 

GRAV  (J  +  NN ) =GR A V (NN-J+1) 

71  H(NN+J)=H(NN)  +  DX*  J 

NN=NN*2 

IF  (JPRINT. LT. 3)  WRITE  (6r509)  (GRAV  (J)  ,  J=1  ,  NN) 

509  FORMAT (/////, 30Xr «  WITH  SEDIMENT  GRAVITY  REMO V  ED  * , //, 1 0 ( 1 
1  *,/)) 

XMAX  =  H  (NN) 

DX=XMAX/NN 

C  OFFSET  EACH  SPLINED  GRAVITY  BY  UM  TO  FACITATE  CONVERGENCE 
DO  7  J=1rNN 
7  GRAV  (J)  =GRAV  (JDM 

IF  (JPRINT.  LT.  3)  WRITE  (6,55)  (H  ( J)  ,  GRA V  ( J)  ,  J=  1  ,  NN) 

55  FORMAT  (/////, 40X, ® DISTA NCE  IN  KM.  AND  SPLINED  GRAVITY  IN  M 
1L*,//, 10 (2F9 .3,6X,  2F9 . 3 , 6X r 2F9. 3 , 6 X , 2F9 . 3, 6X , 2F9 . 3,/) ) 

CC  COMPLETE  THE  INVERSION  OF  THE  GRAVITATIONAL  ANOMALY 
C  INPUT  THE  MAX  NUMBER  OF  TERMS  IN  THE  SERIES  TO  BE  TAKEN. 

C  INPUT  THE  MAX  NUMBER  OF  ITERATIONS  TO  BE  USED. 

C  THE  DESIRED  ACCURACY  BETWEEN  SUCCESSIVE  TOPOGRAPHIC  HERAT 

C  THE  INDEX  IZUP  FOR  ALTERING  THE  VALUE  ZUP 

C  FOR  EACH  ITERATION  IZUP=0  IF  THERE  IS  NO  CHANGE  FOR  ZUP, 

C  ELSE  IZUP  =1.  CHANGING  THE  VALUE  OF  ZUP  (BEING  THE  MEDIAN) 

C  SPEED  UP  THE  ITERATION  PROCESS. 

C  JPRI NT=0  TO  OUTPUT  THE  FOURIER  AMPLITUDE,  JPRINT>2  TO 

C  OUTPUT  THE  LARGEST  ABSOLUTE  VALUE  FOR  THE  FORWARD  ITERATIO 

ITERMS=35 
ITER= 1 5 
AC 1=5. E- 03 
ACC= 1 . E-03 
IZUP=0 

DO  34  J J= 1 , N J J 


CC 

CALCULATE  A  BANDPASS  FILTER 

READ (5,35)  WH,SH 

35 

FORMAT (2F8. 4) 

CALL  LOPASS (B,F, WH,SH,NN,  1) 

WRITE  (6,36)  ( WH, SH ,  (B  (J)  ,J  =  1,NN)  ) 

36 

FORMAT (/// ,20X,23H  LOW  PASS  FILTER 

WH=, F6. 4 , 10X, 4H  SH=, 

1  F  6 . 4  ,  //,  10  (15F8.4,/)  ) 

CC  CALCULATE  THE  LOW  PASSED  GRAVITY  PROFILE 
DO  U2  J=1 , NN 

42  ANOM(J)  =CMPLX  (GRAV  (J)  ,0  .) 
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CALL  FFTTWO(ANOM, NN) 

DO  43  J=1,NN 

43  ANOM(J)  =  ANOM  (J)  *B  (J) 

CALL  FFTTWO (  A  NOM  ,  NN) 

CALL  REORDR  (ANOM,NN) 

DO  44  J=1,NN 

44  DG (J) =REAL (ANOM (J) ) /NN 
WRITE  (6,45)  (DG (J) , J=1, NN) 

45  FORMAT (/////, 30X,31H  FILTERED  GRAVITATIONAL  ANOMALY f //, 1 0 ( 

1  4,/)) 

CC 

cc 

DO  34  KK= 1 , NKK 
ZUP=Z Z  (KK) 

CC 

DO  34  LL=1 ,  NLL 
C  GRHO=6. 6  7*RHO  (LL) 

W RI TE  (6 , 39 )  RHO  (LL)  ,  RHOTWO  (LL) 

39  FORMAT (/////////// //, 20 Xr 1 7H  ASSUMED  DENSITY=/  2F8.2,//) 

DO  40  J=1,NN 

DG4 (J) =0.0 
DG  (J)  =GR AV  (J) 

40  H 1  ( J)  =0.0 
ISKIP=0 
NTERMS=35 

C  CALL  INVERSION  SBROUTI NE  TO  DO  INVERSIO- 

CALL  RHOINV (DG3 r  DG2 , H 1 , H2, H, DG,SUM, HTON, ANOM , B , SH , BTK2, 
1BTK, DERH0,DERH02rDG4,RH0  (LL)  , RHOTWO (LL)  , RODIG1  ,  RODIG2) 
NTERMS=35 

C  REMOVE  OFFSET  AND  ADD  GRAVITY  OF  SEDIMENT 
NN=NN/2 
DO  98  M= 1 , NN 

98  DG(M)  =DG4  (M) +UM-G1  (M) 

CC 

IF  ( J  PRINT.  LT .  3)  WRITE(6,51)  (G2  ( J)  r  DG  ( J)  ,  J=  1 ,  N  N) 

IF  (J PRINT . GE . 3)  WRITE(6#512)  (G2  (J)  , J= 1 r NN) 

512  FORMAT (/////, 30X r '  ORIGINAL  GRA VIT Y 9 r //, 1 0 ( 1 0F  1 2 . 

1  4,/)) 

IF  (JPRINT . GE . 3)  RRITE(6r511)  (DG  ( J)  # J= 1 , NN) 

511  FORMAT (/////» 30X, *  GRAVITY  FROM  INVERTED  TOPOG • , //, 1 0 ( 1  OF 1 

1  *./)) 

51  FORMAT {/////, 40X, 48H  OBSERVED  AND  CALCULATED  GRAVITATIONAL 

1 1 E  S  r // ,  1 0  (2F9. 3 , 6X  r  2F9.3r6Xr  2F9.3f6X,  2F9.3r6X,  2F9. 

34  CONTINUE 

STOP 
END 

SUBROUTINE  RHOBOG (H*DGrSUM,HTON,NTERMS, DELRO, 

1  DEL  R02  r  RODIG2) 

CC$$$$$  CALLS  FFTWO 

C  PROGRAM  COMPUTES  GRAVITY  OF  MODELS  WITH  LATERAL  VARIATION 
C  OF  DENSITY. 

C  DELRO  IS  DENSITY  OF  W  END  OF  MODEL 
C  DELR02  IS  DENSITY  OF  E  END  OF  MODEL 


' 


' 
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C  R0DIG2  STORES  DIGITISED  DENSITY 

CC  COMPUTES  2- D  BOUGUER  ANOMALIES  WITH  FFT  POWER  METHOD 
CC  *  NTERM  S'  IS  THE  NUMBER  OF  TERMS  IN  POWER  SERIES  TO  BE  TAKEN. 
CC  PARAMENTS  COME  THROUGH  COMMON  /HILLS/  *  *NN«  NUMBER  OF  TERMS 
CC  * DX  *  SPACING  IN  X  DOMAIN, «ZUP'  VERTICAL  DISTANCE  ABOVE  Y  =  0  OF 
CC  1 GRHO 1  BIG  G  TIMES  RHO,*H'  ARRAY  WITH  'NN'  ELEMENTS  CONTAININ 
COMMON/HILLS/NN, DX,ZUP, GRHO, GRO 
DIMENSION  RODIG2  (1 )  ,H  (1 ) , DG (1 ) 

COMPLEX  SUM  (1)  ,HTON  (1)  ,  HC 
COMMON/BOUG/ACC, I PRINT 
DATA  TWOPI/6. 2831 85307/ 

DATA  G/6.67/ 

GRO=DELRO*G 

IF  (DELRO. EQ. DELR02)  GOTO  77 
ND=NN/2 

ROIN=( DELRO- DELR02)/(ND-  1) 

DO  303  1=1, ND 

R0DIG2  (I) =DELRO-ROIN*  (1-1) 

RODIG2 (1  +  ND) =DELRO-ROIN*  (ND-I) 

303  CONTINUE 

IF  (IPRINT.LE.  2)  WRITER, 123)  (RODIG2(I)  ,1=1, NN) 

123  FORMAT  (///,' DIGITISED  DE  NSITY  •  , //,  1  0  ( 10  FI  2. 3  , /)  ) 

77  WRITE (6, 55)  DELRO, DELR02 

55  FORMAT (//, 35X, 'DENSITY  =  »,2F8.4) 

C 

WRITE  (6,20)  ACC 

20  FORMAT (///,30X, 29H  ACCURACY  LEVEL  IN  BOUGER  IS = , E 1 2 . 3 , //) 
CC  SETS  VARIOUS  CONSTANTS  AND  ZEROS  ’SUM*  READY  TO  RECEIVE  SERIE 
NNBY2=NN/2 
N2PLUS=NNBY2+ 1 
DK  =  TWOPl/  (DX*NN) 

DO  110  1=1, NN 

110  SUM (I) =(0.0, 0.0) 

F AC= 1 . 0 
BIGM 1=10. 

C 

CC  *  N '  COUNTS  NUMBER  OF  TERMS  IN  SERIES. 

DO  2500  N=1 , NTERMS 
FAC  =  N  *FAC 

CC  COPIES  N-TH  POWER  OF  TOPOGRAPHY  INTO  COMPLEX  ' HTO N  * 

BIG=0 . 0 

IF  (DELRO. EQ. DELP02)  GOTO  211 
DO  233  J= 1 , NN 

2  33  HTON (J) =CMPLX  (RODIG2 ( J) * H ( J) ** N, 0 . 0) 

GOTO  234 

211  DO  210  1=1, NN 

210  HTON  (I)  =CMPLX  (H  (I)  **N  ,  0 . 0) 

& 

234  CALL  FFTT WO  (HTON , NN) 

C 

CC  FILTERS  FOURIER  COEFFS  AND  ADDS  THEM  INTO  'SUM* 

CC  FEWER  THAN  NTERMS  IN  THE  SUM  MAY  BE  TAKEN  IF  AN  ABSOLUTE  CON 
CC  CRITERION  IS  MET. 
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CC  ALSO  PERFORM  THE  UPWARD  CONTINUATION  AT  THIS  STAGE 
TERM=TWOPI*G/NN 

IF  (DELRO. EQ. DELR02)  TERM=TERM*DELRO 
DEXP=EXP  (-DK*ZU  P) 

IF(N.EQ.I)  SUM  (1)  =HTON  (1)  *TERM 

DO  2200  K  =  2 , N2PLUS 

TER M= TERM* DEXP 

C A Y N=  (  (K-1)  *DK) ** (N-  1 ) /FAC 

CAYN=CAYN*TERM 

IF  (K. EQ. N2PLUS)  GO  TO  40 

HC=HTON  (K)  *CA  YN 

B=C ABS (HC) 

IF(B.GT.BIG)  BIG=B 
SUM  (K) =  SUM  (K)  +HC 

40  MI NUS  =  NN-K  +  2 

HC  =  HT  0  N  (  M I N  U  S  ) *CAYN 
B=C ABS (HC) 

IF (B . GT. BIG)  BIG=B 

2200  SUM (MINUS) =SUM  (MINUS)  +  HC 
IF(N.EQ.I)  FI RS/=  BIG 
IF (IPRINT. LE. 5)  WRITE(6,42)  N,  BIG 

42  FORMAT (30X,3H  N  =  ,  1 4  10X,51H  LARGEST  ABSOLUTE  VALUE  OF  TERM 

1WARD  PROBLEM=, E15.4)  . 

BB= AM  AX  1  (BIG, BIGM 1 ) 

I F (  BB/FIRST  .LT.  ACC  )  GO  TO  43 
BIGM 1 =BIG 

2500  CONTINUE 

WRITE  (6,44)  FIRST, BB 

44  FORMAT (////, 5 1H  *****  THE  FORWARD  PROBLEM  HAS  NOT  CONVERGE 
1  , 1  OX, 7H  FIRST=, E  1  2 . 3 , 10X,5H  BIG=,E12.3) 

43  CALL  FFTT WO (SUM , NN ) 

CALL  REORDR (SUM, NN) 

DO  70  1=1, NN 

70  DG  (I)  =REAL  (SUM  (I)  ) 

RETURN 

END 

SUBROUTINE  RHOINV ( DG3 , DG2 , H 1 , H2 , H , DG , SUM, HTON , 

1  A  NOM , F IL, SH, BTK2 , BTK, DER HO, DERH02 , DG4 , RH0 1 , 
1RH02,R0DIG1,R0DIG2) 

C  *####*#****  Jje#**  ***  **  ***  #******:>!£*  *  ****  *  *****  *** 

C  ROUTINE  COMPUTES  THE  TOPOGRAPHY  FROMA  GRAVITATIONAL  ANOMALY 

C  OF  A  MULTILAYER  MODEL,  WITH  LATERAL  VARIATION  OF  DENSITIES 
C  BEDTHICKNESS. 

CC  HI  CONTAINS  THE  INITIAL  GUESSED  TOPOGRAPHY  MEASURED  WITH  RES 

CC  REFERENCE  LINE  ZUP. 

CC  IN  THE  PERTURBATION  EXPANSION  WE  SUM  FROM  2  TO  NTERMS . 

CC  ITER  IS  THE  TOTAL  NUMBER  OF  ITERATIONS  TAKEN  TO  FIND  THE  TOP 
CC  IF  IS KI P=0  THEN  THE  PERTURBATION  TERM  IS  NOT  ADDED. 

CC  IF  THE  MEAN  ROOT  SQUARE  DEVIATION  OF  TWO  SUCCESSIVE  TOPOGRAP 
CC  CALCULATION  IS  LESS  THAN  ACC  THE  PROCEDURE  IS  TERMINATED. 

C  HI  IS  THE  INITIAL  GUESSED  TOPOGRAPHY  WHICH  CAN  BE  INITIALISE 
C  =0.0  IF  THERE  IS  NO  KNOWN  OR  GUESSED  VALUES. 

C  H  IS  A  DUMMY  VARIABLE  FOR  THE  INVERSION  OF  TOPOGRAPHY. 
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C  DG  IS  THE  OBSERVED  GRAVITY. 

C  SOM  IS  THE  SUMMATION  OF  THE  FFT  OF  GRAVITY  AND  TOPO  TERMS. 

C  HTON  IS  A  DUMMY  COMPLEX  VARIABLE  FOR  THE  TOPOGRAPHY  H(X),  WH 

C  FOR  FFT  OF  TOPOGRAPHY. 

C  ANOM  IS  A  DUMMY  VARIABLE  FOR  COMPLEXING  THE  OBSERVED  GRAVITY 
C  FOR  FFT  OF  DRAVITY  DATA. 

C  FIL  IS  THE  FILTER  COEFFICIENTS  SUPPLY  BY  LOPASS  SUBROUTINE. 

C  SH  IS  THE  CUT  OFF  FREQUENCY. 

C  THIS  PROGRAM  EXPECTS  VALUES  FOR  NT  ER  MS  *  ITER, ISKIP, ACCr IZUP, 
C  NNr DX ,ZUP, AND  GRHO  FROM  THE  COMMON  STATEMENT. 

C  IPRINT  IS  THE  INDEX  FOR  PRINTING  THE  VALUES  OF  SUM  FOR  CHECK 

C  ISKIP  IS  THE  INDEX  FOR  SKIPPING  THE  PERTURBATION  TERMS  IF  TH 

C  TOPO  IS  0.  ISKIP  =  0  FOR  TOPO  =0  ELSE  ISKIP  =1 

C  IZUP  IS  THE  INDEX  FOR  CHANGING  THE  VALVE  OF  ZUP  IF  THE  INITI 

C  S  NOT  THE  MEDIAN  OF  THE  TOPO.  UNLESS  THE  ZUP  VALUE  IS  IS  OB 

C  MEDIAN  IZUP  PRETTY  NELL  HAS  TO  BE  =  1 

C  ##  #  #  *  #  #  #  #  *  jjc#  #  #  sje  #  #  #  #  jje##  3{c  ##  #  #  ###  J{c  j{e 

C  GRAVITY  OF  EACH  OVERLYING  LAYER  IS  REMOVED  BEFORE  EACH  ITERA 

C  LINEAR  VARIATION  OF  BEDTHICKNESS  AND  DENSITIES  IS  ASSUMED 

C  IF  DENSITIES  AT  OPPOSITE  ENDS  OF  EACH  LAYER  IS  SAME  DIGITISI 

C  PROCESS  IS  SKIPED.  SAME  IS  TRUE  FOR  BEDTHICKNESS. 

COMMON/HILLS/NN, DX  ,  ZU  P  ,  GRH  0,  G  RO 

COMMON/INV/NTERMS, ITER, ISKIP, ACC ,IZUP,IPRI NT,FRAC, NLAY 

DIMENSION  DG3  (1) ,H2  (1) , DG2  { 1 )  , H 1  { 1 ) , H  (  1 )  , DG  ( 1 )  r FI L  ( 1 ) 
DIMENSION  RODIG1  (1)  ,RODIG2 (1)  , DERH02  (1)  f  BTK  (1) , BTK2  (1)  ,  DERH 
DIMENSION  DG4(1) 

Q  ###******#****#*#**#*#*&#******#****************** 

COMPLEX  SU  M  (  1 )  ,  HTON  (1)  , A  NOM ( 1 ) 

COMPLEX  HC 

DATA  PI/3.  lui5926535/,G/6.67/ 

CC 

TWOPI=2*PI 
NND2=NN/2 
N2PLUS=NND2+1 
DK  =  TW  PPI/  (DX*NN) 

CC 

CC  CALCULATE  CONTEIBUTIONS  BELOW  THE  CUTOFF  FREQUENCY 
NFINAL=NN*SH  +2 

IF (NFINAL  -GT.  N2PLUS)  NFI NAL=N2 PLUS 
STORE  INPUT  GRAVITY  IN  DG2  TO  FIND  RESIDUE  GRAVITY 
A  FRACTION  OF  INPUT  GRAVITY  IS  USED  FOR  FIRST  ITERATION 
DG3  STORE  GRAVITY  TO  FIND  THE  AVERAGE  OF  TWO  ADJ  ITERATIONS 
DO  9  1=1, NN 
DG2  (I)  =DG  (I) 

DG  (I) =DG  (I)  *FRAC 
9  DG3  (I)  =DG  (I) 

IF  (IPRINT. LE. 2)  WRITE  (6,8)  (DG (I) , 1= 1 , NN) 

8  FORMAT (//,20X, 'GRAVITY  USED  FOR  FIRST  INVERSION', 

1//,  10(10F12.3,/)  ) 

DO  10  J=1,NN 
SUM (J) =CMPLX  (0. , 


0.) 
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10  ANOM(J)  =CMPLX  (DG(J)  ,0.) 

6  FORMAT (//,  1  OFFSET  GRAVITY  FOR  IN VERSION  * , FI  0 . 5, //, 

1 10 (10F12. 3,/) ) 

C 

CALL  FFTTWO  (ANOM,  NN) 

C 

SOL  D= 1 . E0  5 
C 

DO  50  JJ=1 , ITER 
C 

CC  CONVERGENCE  IS  MOST  RAPID  IF  THE  TOPOGRAPHY  IS  EQUALLY  POSIT 
CC  WITH  RESPECT  TO  ZUP.  FIND  THE  LARGERS  AND  SMALLEST  VALUES  0 
C 

BIG=-1. E09 
SMALL= 1 . E09 
DO  12  J=1,NN 

IF  (HI  (J)  .  GT.  BIG)  BIG=H1  (J) 

IF  (HI  (J)  .LT.  SMALL)  SMALL=H1(J) 

12  CONTINUE 

DELZU  P=- (BIG+SMALL)  /2 
C  COMPUTE  ELEVATION  OF  TOP  HORIZON 

WRITE  (6, 1 3)  JJ, BIG, SMALL, ZUP, DELZU P 

13  FORMAT (//, 5X, 6H  ITER= , 14, 5X, 1  OH  MAX  EL E V= , F 1 0 .  4 , 5X , 1  OH  MIN 

1  F10.4,5X,5H  ZUP=,F8. 4,5X,8H  DELZUP=, F8. 4) 

IF  (IZUP.EQ.O)  GO  TO  31 

ZUP=ZUP+DBLZUP 

DO  14  J=1,NN 

14  HI  (J) =H1  (J) +DELZUP 
C 

CC  INITIALIZE  SUM 

Q  *  ##***:<£****  **:{£***  *********** 

31  DEXP=EXP (DK*ZUP) 

TERM=1./(TWOPI*G) 

C  IF  DENSITIES  AT  OPPOSITE  ENDS  ARE  SAME  USE  CONSTANT  DENSITY 
IF  ( RH01 . EQ. RH02)  TERM=TERM/RH01 
SUM  (1  )  =  ANOM  (1)  *TERM 
DO  15  K=2, NFINAL 
TER M= TERM* DEXP 

IF  (K. NE. N2PLUS)  SUM(K)  =  ANOM  (K)  *TERM 
MI NUS=NN-K+2 

15  SUM  (MINUS) =  ANOM  (MINUS)  *TE RM 
****:{£  #**:{::*:{£************  ****** 

IF  (IPRINT. EQ. 0)  WRITE  (6,55)  (SUM  ( J)  ,  J  = 1 , NN) 

55  FORMAT (////, 19H  FOURIER  AMPLITUDES,//,  1 0  ( 1 0 E 1 2 . 3 , /) ) 

C  DIGITISE  DENSITY 

IF  (RHO 1 . EQ. RH02)  GOTO  621 
ND2=  NN/2 

ROIN= (RH01-RH02) / (ND2- 1) 

DO  303  1=1, ND2 

RODIG1  (I) =RH01-R0IN* (1-1 ) 

R0DIG1  (I  +  ND2) =  RH0 1 - ROI N*  (ND2-I) 

303  CONTINUE 


. 
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C  ******************  **************** 

IF  (IPRINT.LE.  2)  WRITE(6,373)  (RODIG 1  (I )  ,1=  1 ,  NN) 

373  FORMAT (///,10Xr •  DIGITISED  DE NSIT Y= 1  , 

1  //,  1  0  ( 10F 1 2 . 4 , /)  ) 

C  IF  THE  INITIAL  TOPOGRAPHU  IS  ZERO  THEN  DO  NOT  CALCULATE  THE 
621  IF  (ISKIP. EQ. 0)  GO  TO  20 
FAC=1 .0 

DO  16  N=2 , NTERMS 
FAC=N*FAC 

CC  COPY  THE  NTH  POWER  OF  THE  TOPOGRAPHY  INTO  COMPLEX  HTON 
IF  ( RHOI . EQ. RH02)  GOTO  217 
DO  724  J=1,NN 

7  24  HTON  (J) =CMPLX (ROD TGI  (J) *H1  (J) **N, 0. 0) 

GOTO  127 

217  DO  17  J  =  1 r  N  N 

17  HTON(J)  =CMPLX  (HI  (J)  **N,  0.) 

C 

127  CALL  FFTTWO (HTON,NN) 

C 

c 

CC  FIND  THE  LARGEST  VALUE  OF  THE  PERTURBATION  TERM 
BIG=0 .0 

DO  18  K=2  , NFINAL 

CAYN=  l  (K- 1 J *DK) ** IN- 1)  /FAC 

IF  (K.  EQ.  N2PLUS)  GO  TO  40 

HC=  HTON (K)  *CAYN 

B=C ABS  (HC ) 

IF  ( B . GT . BIG)  BIG  =  B 
SUM  (K)  =SUM  (K-HC 

40  MINUS  =  NN-K+  2 
HC=HTON (MINUS) *CAYN 
B=C ABS (HC) 

IF (B. GT. BIG)  BIG=B 

18  SUM  (MINUS)  =SUM  (MINUS)  -  HC 
IF(N.EQ.2)  FI RST=  BIG 

IF  (IPRINT.LE. 2)  WRITE(6,41)  N,BIG 

41  FORMAT (30X,3H  N =  , 1 4 , 1  OX r 5 1 H  LARGEST  ABSOLUTE  VALUE  OF  TERM 
1WARD  PROBLEM^ ,  E15.4) 

I F ( B I G /F I R S T  .LT.  5. E-03)  GO  TO  20 
16  CONTINUE 

WRITE  (6,56) 

56  FORMA T(20X,34H  THE  FORWARD  SUM  HAS  NOT  CONVERGED) 

C 

CC  APPLY  THE  LOWPASS  FILTER 

20  IF  (IPRINT.EQ.  0)  WRITE(6,55)  ( SUM  ( J)  , J= 1 , NN) 

IF  (RHOI . EQ. RH02)  GOTO  521 

C  ****************************** 

C  INVERSE  TRANSFORM  TO  GET  TOPOG,  DIVIDE  BY  CORRESPONDING 

C  DENSITY,  TRANSFORM  AND  FILTER  AND  COMPARE  EACH  ADJACENT 

C  ITERATION. 

CALL  FFTTWO (SUM,NN) 

CALL  REORDR (SUM , NN) 

DO  272  J= 1 , NN 


' 
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2  72  H (J)  =RE AL  (SUM (J)  ) /  (NN*  RODIG 1  (J) ) 

DO  299  J=1,NN 

299  SUM(J)  =CMPLX  (H  (J)  ,0.0) 

CALL  FFTTHO (SUM,NN) 

DO  37  K=1,NN 

37  SUM  (K)  =SUM  (K)  *FIL(K) 

CALL  FFTTWO (SUM,NN) 

CALL  REORDR (SUM,  NN) 

S1=0.0 

DO  263  1=1, NN 
H  (I)  =REAL  (SUM  (I)  )  /NN 
H2 (I) =H  (I) 

263  S  1  =  S  1  + ( H (I- H 1 (I)  )  **  2 
C  **#*###  #**#*##*:  #######*## 

GOTO  727 

521  DO  32  J=1,NN 

32  SUM  (J)  =SUM  (J)  *FIL  (J) 

IF (IPRINT. EQ. 0)  WRITE  (6,55)  (SUM  ( J)  , J= 1 , NN) 

C 

CC  CALCULATE  THE  TOPOGRAPHY  BY  TAKING  THE  FOURER  TRANSFORM 
CALL  FFTTWO (SUM, NN) 

CALL  REORDR  (SUM, NN) 

C 

S1=0.0 

DO  22  J= 1 , NN 
H  (J)  =  RE  AL  (SUM  (J)  )  /NN 
H2 (J) =H  (J) 

22  S1  =  S1+ (H  (J-H1  (J)  )  **2 

C  *#***##  #****#  Jit*#*##*##;!}:# 

727  S1  =  SQRT  (S1)/NN 

WRITE (6,23)  JJ,S1 

23  FORMAT (//, 40 X , 1 1 H  ITE RA TION= , 15 , 1  OX , 7 H  ERROR  =  , E12. 4) 

IF  (IPRINT.  LE.  2)  WRITE(6,24)  ZUP,  (H  (  J)  ,  J=  1  ,  NN) 

24  FORMAT (///, 10X,41H  THE  TOPOGRAPHY  MEASURED  RELATIVE  TO  ZUP 

1  F8.4,//,  10(1  OF  12.  4,/)  ) 

C 

CC  CHECK  TO  SEE  IF  ITERATION  IS  STARTING  TO  DIVERGE 
’  XX=  (SI-SOLD)  /SOLD 
C 

C  DO  NOT  CHECK  FOR  DIVERGENCE  OF  SECOND  ITERATION  BECAUSE 

C  FIRST  ITERTN  USE  VERY  DIFFERNT  GRAVITY 

IF  (JJ.EQ.2)  GOTO  42 
IF  (XX  .  LT.  .01)  GO  TO  42 
WRITE  (6,43) 

43  FORMAT (///,20X,57H  ****  THE  ITERATION  PROCEDURE  IS  STARTIN 
1  VERGE  ****,///) 

USE  THE  TOPOGRAPHY  FROM  THE  LAST  ITERATION 
DO  44  J= 1 , NN 
H2  (J  )  =H  1  (J) 

44  H  ( J)  =H1  (J) 

GO  TO  30 

C  ****#********#*:*❖**********************  ********** 
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42  IF (SI .LE. ACC)  GO  TO  30 
IF  (JJ.EQ.ITER)  GOTO  381 
SOLD=S1 

********************  j^******.^,^^ 

IF  ANY  OVERLYING  LAYER  HAS  ZERO  BEDTHICKNESS  (I. E. )  TWO  LAYER 
DO  NOT  COMPUTE  GRAVITY  OF  OVERLYING  LAYERS. 

DO  711  K=1,NN 
711  DGU(K)=0.0 

DZUP=ZUP 
N  ND2  =  N  N/2 
NBD=NL AY-2 

C  IF  BED  IS  OF  SAME  THICKNESS  DO  NOT  CALCULATE  VARIATION 
DO  909  1= 1 , NBD 

IF  (BTK(I)  ,EQ. BTK2  (I)  . AND.BTK(I)  .EQ.0.0)  GOTO  742 
IF  (BTK(I)  .EQ.BTK2  (I)  )  GOTO  623 
DBTK=BTK (I-BTK2 (I) 

BTKIN=DBTK/ (NND2-1) 

DO  99  L= 1 r  NND2 

H2  (L)  =H2  (L-BTKIN*  (L-1) 

9  9  H2  ( NND  2  +  L)  =H2  (NND2+L-BTKIN*  (NND2-L) 

C  READJUST  HORZN  RESPONSIBLE  FOR  CORRESPONDING  BED 
623  ZUP=ZUP-BTK (I) 

IF  (IPRINT.LE.  2)  WRITER, 75)  Z  UP  ,  (H2  ( J )  ,  J=  1  ,  NN) 

75  FORMAT (///, 1 0  Xr  4 1 H  THE  TOPOGRAPHY  MEASURED  RELATIVE  TO  ZUP 

1  F8.  4r //,  10  (10F12. 4,/) ) 

C  COMPUTE  GRAVITY  OF  OVERLYING  LAYER 

CALL  RHOBOG(H2,DG,SUM,HTON,NTERMS,DERHO(I)  ,DERH02(I)  , RODIG2 
IF  (IPRINT.LE. 2)  WRITE(6,74)  ZUP  ,  (DG  (J)  ,  J=  1  ,  NN ) 

DO  919  K=1,NN 
919  DG4  (K)  =DG4  (K)  +DG  (K) 

909  CONTINUE 

q  ***************************************** ******** 

74  F 0 R M A T ( /// ,  1 0  X  ,  *  THE  GRAVITY  CALCULATED  RELATIVE  TO  ZUP=«, 
1  F8.4,//,  10  (10F12. 0,/) ) 

C  USING  THE  AVG  OF  THE  GRAVITY  USED  IN  THE  LAST  ITERATION  AN 

C  THE  DIFFERENCE  OF  THE  INPUT  GRAVITY  AND  THE  GRAVITY  OF 

C  OVERLYING  LAYER 

C  COMPUTE  THE  RESIDUAL  GRAVITY 

DO  88  M=1,NN 

DG  (M)  =  (DG2  (M-DG4  (M)  +DG3  (M)  )  /2 . 0 
88  DG3  (M)  =DG  (M) 

IF  (IPRINT.LE. 2)  WRITE(6,76)  DZ  UP ,  (DG  (J )  ,  J=  1  r  N  N) 

76  FORMAT (///, 10Xr 1  THE  RESIDUE  GRAVITY  CALCULATED  RELATIVE  T 

1  ,F8.4r//,  10  (10F12.4,/)  ) 

C  RESTORE  ZUP  TO  ITS  STARTING  VALUE 
ZUP=DZUP 
DO  77  J= 1 r  N  N 
SUM  (J)  =CMPLX  (0.  ,  0. ) 

77  ANOM  (J)  =CMPLX  (DG  (J)  ,  0.  ) 

CALL  FFTTWO  (ANOM,NN) 

CC  TRANSFER  THE  TOPOGRAPHY  BACK  TO  HI 
742  DO  25  J=1,NN 
25  HI  (J) =H  (J) 


' 
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50  I SK IP= 1 

C 

381  WRITE  (6,26)  SI,  ACC  ,  ITER 

26  FORMAT (///,44H  *********  CONVERGENCE  HAS  NOT  BEEN  OBTAINED 
1  10Xr7H  ERROR=,E12. 3,1QX,18H  DESIRED  ACCURACY= , E 1 2 . 3 , 1  OX , 

1  22H  NUMBER  OF  ITERATIONS^, 15) 

C  COMPUTE  FINAL  INVERTED  GRAVITY  OF  INVERTED  TOPOG. 

30  MN=NN/2 
ZUP=DZUP 

IF  (IPRINT.GT. 2)  WRITE  (6,24)  ZUP,  (H  (J)  , J=1 , MN) 

CALL  RHOBOG(H2rDG,SUM, HTON, NTERMS , RHOI , RH02 , RODIG2 ) 

DO  633  N=1rMN 
DG4  (N)  =0.0 
633  DG4  IN)  =DG  (N) 

IF  (IPRINT.LT. 2)  WRITE  (6,7  4)  ZUP ,  (DG  ( J )  , J= 1 , MN ) 

DO  809  1=1, NBD 

IF  ANY  OVERLYING  LAYER  HAS  ZERO  BEDTHI CK NESS (I . E. )  T WO  LAYER 
DO  NOT  COMPUTE  GRAVITY  OF  OVERLYING  LAYERS. 

IF  (BTK (I)  .EQ. BTK2 (I)  . AND.BTK  (I)  . EQ.0.0)  GOTO  809 
IF  (BTK (I)  . EQ. BTK2  (I) )  GOTO  624 
DBTK=BTK (I-BTK2 (I) 

BTKI N=  DBTK/ (NN  D2- 1 ) 

DO  701  L= 1 , NN D2 
H2  (L)  =H2  (L-BTKIN*  (L - 1 ) 

701  H2  (NND2+L) =H2 (NND2+L-BTKIN* (NND2-L) 

624  ZUP=ZUP-BTK  (I) 

IF  (IPRINT.GT.  2)  WRITE(6,24)  Z  UP  ,  (H2  (J)  ,  J=  1  ,  MN) 

CALL  RHOBOG (H2,DG, SUM, HTON, NTERMS, DERHO(I)  ,DERH02(I)  ,RODIG2 
IF  (IPRINT.LT. 2)  WRITE(6,74)  ZUP ,  (DG  (J)  , J=  1  , MN) 

DO  929  K= 1 , MN 
929  DG4  (K)  =DG(K)  +DG4  (K) 

809  CONTINUE 
RETURN 
END 

SUBROUTINE  LOPASS  (B, F, W H,SH, N, II) 

CALCULATE  A  HIGH  PASS  FILTER  USING  A  COSINE  TAPER  BETWEEN  WH 
AT  FREQUENCIES  LESS  THAN  WH  THE  FILTER  IS  EQUAL  TO  UNITY,  AND 
GREATER  THAN  SH  THE  FILTER  IS  ZERO. 

IF  11=1  A  FULL  TWO  SIDED  FILTER  COMPATIVLE  WITH  FOURERR  TRANS 

DIMENSION  B  (1)  ,F  (1) 

DATA  PI/3.1415926535/ 

NN=N/2+ 1 
DO  11  J= 1 , NN 

1.4  IS  MULTIPLIED  BECAUSE  DIFF  SAMPLING  INTERVALS  WERE  USED 
5  KM.  AND  7  KM.  WERE  USED  IN  TWO  SEPARATE  PROFILES. 

11  F(J)  =  ( J-  1 )  *1.4/(N*1.0) 

DO  10  J= 1 , NN 
IF  (F  (J-WH)  20,20,21 
20  B (J)  =1 .0 

GO  TO  10 

IF  (F  (J -SH)  24,24,25 


21 


' 
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24  B  (J)  =  (1.  -  COS  (PI*  (F  (J-SH) /(WH-SH)  ))  /2  . 

GO  TO  10 

25  B  (J)  =0.0 

10  CONTINUE 

C 

IF  (II. NE.  1)  RETURN 

C  CALCULATE  A  TWO  SIDED  FILTER  SERIES 
NN=N/2-1 
DO  30  J=1,NN 
30  B  ( N+  1-J)  =  B  (J-s- 1 ) 

RETURN 

END 

SUBROUTINE  SPLNCO  (N1,N2rXrY,B,C,D) 

C  COMPUTATION  OF  THE  COEFFICIENTS  OF  A  CUBIC  SPLINE 

C  INTERPOLATING  BETWEEN  GIVEN  DATA  POINTS. 

C  INPUT. 

C  N 1 , N2  NUMBER  OF  FIRST  AND  LAST  DATA  POINT (0 . LT . N1 . LT. N2) 

C  X(I)  rY  (I)  ,I  =  N1  ,N1  +  1  ,.  .  .  ,N2  ARRAYS  WITH  X(1)  AND  Y(1)  AS  ABS 

C  AND  ORDINATE  OF  I-TH  DATA  POINT.  THE  COMPONENTS  OF  THE 

C  ARRAY  X  MUST  BE  EITHER  STRICTLY  MONOTONIC  INCREASING 

C  OR  DECREASING. 

C  OUTPUT. 

C  B (I)  ,C  ll)  , D  (I)  rI=N1,N1  +  1,..,N2  ARRAYS  COLLECTING  THE  COEFFI 

C  THE  CUBIC  SPLINE  F(XX).  IF  XX  LIES  BETWEEN  Xll)  AND  X( 

C  THEN  F  (XX)=  (  (D  (I)  *H+C(I)  )  *H  +  B  (I)  )  *H«-Y  (I)  , 

C  WHERE  H=XX-X  (I) . 

C  FURTHERMORE ,C  (N2) =0  WHILE  B(N2)  AND  D(N2)  ARE  LEFT 

C  UNDEFINED. 

DIMENSION  X(1),Y(1),B(1),C(1),D(1) 

M1=N1+  1 
M2  =  N2-  1 
S=0. 

M3  =  M  1  +  M  2 
DO  1  K  =  N 1 , M  2 
D (K)  =X (K+ 1-X  (K) 

R=  (Y  (K+  1  -Y  (K)  )  /D  (K) 

C  (K)  =R~S 

1  S  =  R 

C  (N2) =0 
C  (N1) =0 
S  =  0 
R=  0 

DO  2  K=M1,M2 
C (K)  =C (K) +  R*C  i K~ 1 ) 

B  (K)  =lX(K-1-X(K+1)  )  *2 .  -R  *S 
S=D  (K) 

2  R=S/B(K) 

DO  3  K  =  M  1 , M2 
L=M3-K 

C  (L)  =(D  (L)  *C  (L  +  1-C  (L)  )  /B  (L) 

DO  4  K=N1,M2 

B  (K)  =  (Y  (K+1-Y  (K)  )  /D  (K-  (C  (K)  +C  (K)  +C  (K+1)  )  *D  (K) 

D  (K)  =  (C  ( K+ 1  -C  (K)  )  /D  (K) 


3 


in  \o 


131 


4  C(K)=3.*C(K) 

RETURN 

END 

SUBROUTINE  SPLNEV  ( N 1 , N2 , X , Y , B , C , Dr XI NT , YINT , I) 

C  EVALUATION  OF  A  CUBIC  SPLINE  AT  ABSCISSA  XINT. 

C  INPUT  DATA. 

C  N1,N2  NUMBER  OF  FIRST  AND  LAST  DATA  POINT  lO.LT.N1.LT.N2 

C  X  (K)  /  Y(K)  , K=N1,N1+1,N2  ARRAYS  WITH  X  (K)  AND  Y  (K)  AS  ABSCISS 

C  AND  ORDINATE  OF  THE  K-TH  DATA  POINT.  THE  COMPONENTS 

C  OF  X  MUST  BE  EITHER  STRICTLY  M ONQTONIC  INCREASING  OR 

C  DECREASING. 

C  B  (K)  ,  C  (K)  ,  D  (K)  ,K=N1,N1+1,..,N2  SPECIFY  THE  SPLINE  FUNCTION 

C  TO  BE  EVALUATED.  IF  XX  LIES  BETWEEN  X(K)  AND  X(K+1), 

C  (I.E  .  WITHIN  THE  K-TH  INTERVAL) ,THEN 

C  F  (XX)  =  l  (D  (I)  *H  +  C  (I)  )  *H+B  (I)  )  *H+Y  (I) 

C  WHERE  H=XX- X (I)  .  . 

C  XINT  ABSCISSA,  FOR  WHICH  YINT=F (XINT)  HAS  TO  BE  EVALUATED. 

C  I  IF  ON  ENTRY  1=0,  THEN  THE  NUMBER  OF  THE  INTERVAL,  IN  WHI 

C  XINT  LIES,  IS  NOT  KNOWN  (ERROR  MESSAGE,  IF  XINT  LIES 

C  NOT  BETWEEN  X(N1)  AND  X  (N2)  )  . 

C  ON  LEAVING  THE  ROUTINE,  I  IS  EQUAL  TO  THIS  INTERVAL 

C  NUMBER. 

C  IF  ON  ENTRY  I.NE.O,  THEN  I  SPECIFIES  THE  NUMBER  OF  THE 

C  INTERVAL,  WHERE  XINT  IS  LOCATED  (ERROR  MESSAGE  IF 

C  I.LT.O,  0.LT.I.LT.N1  OR  N2.LT.  I.  UPON  EXIT  ,  I  IS  NOT 

C  OUTPUT  DATA. 

C  I  (SEE  ABOVE) 

C  YINT  =  F  (XINT)  . 

C 

DIMENSION  X  (1)  ,Y  (1)  ,B  (1)  ,C  (  1)  ,  D  (1) 

IF  (I)  16,1,11 

11  IF(I.LT.N1  .  OR.  I.  GE  **N2)  GO  TO  16 
GO  TO  12 

1  IS=1 

IF  (X  (N  1  - X  (N 2)  )  2,3,3 

3  IS=-IS 

2  IF(  (XINT-X  (N1)  )*IS)  10, 14, <i 

14  YINT=Y  (N 1 ) 

RETURN 

4  IF ( (X (N2-XINT) *IS)  10,15,9 

15  YINT=Y(N2) 

RETURN 

9  I  =  N  1 

I2  =  N2 

IF  (I2-I-1)  12,1 2,6 
13= (I+I2+ . 1) *. 5 
IF  (  (XINT-X  (13)  )*IS)  7,13,8 

7  12=13 
GO  TO  5 

8  1=13 

GO  TO  5 

13  YINT=Y  (13) 

RETURN 
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12  H=XINT-X  (I) 

YINT=  l  ID  (I)  *H+C(I)  )  *H+B  (I)  )  *H  +  Y  ( I ) 

RETURN 

10  PRINT  1 0 1 , XINT 

101  FORMAT  (8H  XINT  =  E18.9,32H,  OUT  OF  RANGE  FOR  INTERPOLATION) 
RETURN 

16  PRINT  102,1 

102  FORMAT (5H  I  =  I5,32H,  OUT  OF  RANGE  FOR  INTERPOLATION) 

RETURN 

END 

SUBROUTINE  FFTTWO  (A , N) 

C 

CC  ROUTINE  FOR  FINDING  DISCRETE  FOURIER  TRANSFORM  OF  A  COMPLEX 
CC  MODIFIED  FROM  COOLEY  ET  AL  IEEE  TRANS  OF  EDUCATION  1 2  NO  1  ( 

CC  ARGUMENT  LIST 

CC  A  IS  A  COMPLEX  ARRAY  CONTAINING  THE  SERIES  TO  BE  TRANSFORMED 

CC  RESULTS  ARE  RETURNED  IN  A  0 VER-WRITEING  THE  ORIGINAL  SERIES 

CC  N  IS  THE  NUMBER  OF  TERMS  IN  THE  SERIES.  IT  MUST  BE  A  POWER 

CC  THE  ROUTINE  WILL  PRINT  A  MESSAGE  AND  RETURN  WITHOUT  DOING  AN 

CC  FINDS  FOR  M=0,N-1  SUM,  A  (K  +  1)  *EXP(2*PI*  (0,  1)*  (K+1)  *M+1)  )  ,  WIT 

CC  K  =  0 , N- 1  . 

CC  NOTE  INDIECES  ARE  DISPLACED  WITH  RESPECT  OF  FREQUENCY,  EG.  Z 
CC  FREQUENCY  IS  HELD  IN  A(1). 

C 

DIMENSION  A  (N) 

COMPLEX  A  ,  U  ,  W  ,T 
DATA  PI/3.1415926535/ 

M=ALOG  (FLOAT  (N)  )/0. 693147 
IF (N. NE. 2**M)  GO  TO  1000 
NV  2=N/2 
NM 1=N- 1 
J=1 

DO  7  1=1, NM1 
I F ( I. GE. J)  GO  TO  5 
T= A  ( J ) 

A  (J  )  =  A  (I) 

A  (I)  =T 

5  K=NV2 

6  I F (K . GE. J)  GO  TO  7 

J= J  -K 

K=K/2 
GO  TO  6 

7  J=J+K 
LE=  1 

DO  20  L= 1 , M 
LEI =LE 
LE=2*LE 
U= (1. ,0. ) 

B=PI/FLO AT (LEI) 

W=C MPLX (COS (B) ,SIN (B)  ) 

DO  20  J= 1 , LEI 
DO  10  1= J , N, LE 
IP=I+LE 1 


• 
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T  =  A  (IP)  *U 
A  (IP)  =  A  (I-T 
10  A  (I)  =  A  (I)  +T 

20  U=U*W 

RETORN 

1000  PRINT  1  00  ,  N 

100  FORMAT (5  (4H  ***),I6,70H  IS  NOT  A  POWER  OF  2.  FFTWO  CANNO 
1  FORM  A  SERIES  OF  THIS  LENGTH  ) 

RETURN 

END 

SUBROUTINE  REOR DR  ( A  ,  N) 

COMPLEX  A  (1)  f  DUM 

CC  ROUTINE  REORDERS  THE  OUTPUT  FROM  A  FOURERR  TRANSFORM  ROUTINE 
ND  2=N/2 
DO  10  J=  2, ND2 
DUM  =  A  (J) 

A  (J) = A  (N  +  2-J) 

10  A  (N  +  2- J) =DUM 

RETURN 
END 


' 
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Program  to  invert  residual  gravity  to  obtain  density. 


C  THIS  PROGRAM  INVERTS  RESIDUAL  GRAVITY  TO  OBTAIN  A  DENSITY 

C  -CONTRAST  PROFILE. 

C  WRITEN  BY  JOSEPH  S.K.  LEE,  UNIVERSITY  OF  ALBERTA,  1977. 

C  THIS  PROGRAM  CALLS  SUBROUTINE  DENSTY  AND  DENBOU,  AND 

C  ASSUMES  INPUT  RESIDUAL  GRAVITY  TO  BE  EQUALLY  SPACED. 

C  IT  IS  TO  BE  NOTED  THAT  THIS  PROGRAM  WAS  MODIFIED  FROM  THE 

C  MULTILAYER  INVERSION  SCHEME,  AND  NO  TIME  WAS  SPENT  TO  CHANG 

C  THE  VARIABLES  TO  MEANINGFUL  VARIABLES. 

C  INPUT  TO  THE  MAIN  PROGRAM  SHOULD  BE  OF  THE  FOLLOWING  FORMAT 

C  **  NJJ, NLL, NKK,NN,IG, ITOP, JPRINT.  (#  OF  FILTERS,  #  OF  SLABS, 

C  #  OF  HORIZONS  Z,  TOTAL  NUMBER  OF  RESIDUAL  GRAVITY  STATIONS 

C  FOR  GAVITY  IN  TENTH  OF  MGAL  IG  .NE.  1,IF  H (X)  ALONG  PROFIL 

C  IS  GIVEN  ITOP.EQ.  1  ELSE  ITOP  .NE.  1 

C  JPRINT  .EQ.  3  FOR  NO  OUTPUT  OF  INTERMEDIATE  RESULTS. 

C  **  DX  (SAMPLING  INTERVAL  OF  RESIDUAL  GRAVITY) 

C  **  RHO  (THKNES  OF  SLAB  AT  W  END  OF  PROFL  FOR  ITOP . N E. 1 ; NLL  TERM 
C  **  RHOTWO  (SIMILAR  TO  RHO,  BUT  AT  THE  E  END;  NLL  TERMS) 

C  **  RODIG1  (DIGITISED  H (X)  ALONG  PROFL, FOR  ITOP.EQ.  1,  NN  TERMS) 
C  **  GG  (RESIDUAL  GRAVITY,  NN  TERMS) 

C  **  WH, SH  (FILTER  PARAMETERS,  NJJ  TERMS) 

C  *  ************  *  ************************************** 

DIMENSION  HI  (512)  ,H  (512)  ,DG  (512)  ,  RHO  (512)  ,  RHOTWO  (512) 
DIMENSION  RODIG1  (512)  ,F(512)  ,ZZ  (512)  ,B  (512) 

DIMENSION  A A  (512) ,BB  (512) ,CC (512)  ,GG  (512) 

COMPLEX  SUM  (512),  HTON  (512)  ,  ANOM  (512) 

COMMON /BO UG/ ACC, I PRINT 
COMMON/HILLS/NN, DX , ZUP , GRHO, GRO 

COMMON/IN V/ITERMS, ITER, ISKIP, AC1 , IZUP, JPRINT, FR AC, NLA Y 
DATA  PI/3. 1 415926535/ 

DATA  G/6.67/ 

READ (5, 100)  NJJ,NLL,NKK,NN, IG, ITOP, JPRINT 
READ (5,171)  DX 

100  FORM  AT  (1014) 

101  FORMAT  (10F10. 3) 

171  FORMAT (2F8. 3) 

READ  (5,  10  2)  (ZZ(I)  ,I=1,NKK) 

IF  (ITOP  .NE.  1)  READ  (5, 102)  (RHO (I)  ,1  =  1 , NLL) 

IF  (ITOP  .NE.  1)  READ  (5,  102)  ( RHOTWO  (I)  ,1=  1  ,  NLL) 

IF  (ITOP  .EQ.  1)  READ  (5,47)  (RODIG 1  (J)  ,  J=  1  ,  NN) 

IPRINT= JPRINT 

102  FORMAT  (5F8. 3) 

C  INPUT  THE  #  OF  TERMS  TO  BE  USED  IN  POWER  SERIES 
NTERMS=35 

47  FORMAT (10F12. 3) 

READ  (5,4  7)  (GG  ( J)  ,  J=  1  ,  NN) 

NBY2=NN 

NN=NN*2 


. 
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C  GENERATE  MIRROR  IMAGE. 

DO  37  1=1 , NBY2 
37  GG  (NN-I+1)  =  GG(I) 

IF (ITO P  .NE.  1)  GO  TO  55 
DO  49  J=  1  ,  NBY2 

ti9  RODIG1  (NN-J+1)  =RODIG1  (J) 

C  IF  INPUT  GRAVITY  IS  TENTH  MILLIGAL  CONVERT  TO  MILLIGAL 

55  IF  (IG.EQ.1)  GO  TO  99 
DO  1  1=1 ,NN 

1  GG(I)  =GG  (I)  #0.1 

11  FORMAT  ( 10F1  2.  1) 

99  IF  (JPRINT  .LT.  3)  WRITE  *6,56)  (GG  (J)  , J=  1  r NN) 

IF  (ITOP  .EQ.  1)  WRITE  (6,65)  (RODIG  1  ( J)  r  J=  1  ,  NN) 

65  FORMAT  (/////, 40X,'  INITIAL  TOPOG*, 

1  //,  10 (10F9.3,/) ) 

56  FORMAT  (/////, 40X,*  INITIAL  GRAVITY  ANOMALY*, 

1  //, 10 ( 10F9. 3,/) ) 

CC  COMPLETE  THE  INVERSION  OF  THE  GRAVITATIONAL  ANOMALY 
C  INPUT  THE  MAX  NUMBER  OF  TERMS  IN  THE  SERIES  TO  BE  TAKEN. 

C  INPUT  THE  MAX  NUMBER  OF  ITERATIONS  TO  BE  USED. 

C  THE  DESIRED  ACCURACY  BETWEEN  SUCCESSIVE  TOPOGRAPHIC  ITERAT 

C  DISREGARD  STAMENTS  CONCERNING  IZUP  AND  ZUP,  THEY  HAVE  NO  CONC 
C  THIS  PROGRAM. 

C  THE  INDEX  IZUP  FOR  ALTERING  THE  VALUE  ZUP 

C  FOR  EACH  ITERATION  IZUP=0  IF  THERE  IS  NO  CHANGE  FOR  ZUP, 

C  ELSE  IZUP  =1.  CHANGING  THE  VALUE  OF  ZUP  (BEING  THE  MEDIAN) 

C  SPEED  UP  THE  ITERATION  PROCESS. 

C  J PRI NT  =  0  TO  OUTPUT  THE  FOURIER  AMPLITUDE,  JPRINT>2  TO 

C  OUTPUT  THE  LARGEST  ABSOLUTE  VALUE  FOR  THE  FORWARD  ITERATIO 

ITERMS=35 
ITE  R  =  1 5 
AC 1=5. E-05 
ACC= 1 . E-0  3 
DO  34  J J= 1 , N J J 

CC  CALCULATE  A  BANDPASS  FILTER 
READ(5,35)  WH,SH 

35  FORMAT (2F8. 4) 

CALL  LOPASS  (B,F,WH,SH,NN, 1) 

IF  (JPRINT  .LT.  3)  WRITE(6,36)  WH , SH,  (B  ( J)  , J= 1 , NN) 

36  FORMAT (///,20X,23H  LOW  PASS  FILTER  W H= , F6 . 4 , 1 0 X , 4 H  SH=, 

1  F 6 . 4  , // ,  10  (15F8.4,/)  ) 

CC  CALCULATE  THE  LOW  PASSED  GRAVITY  PROFILE 
DO  42  J= 1 , NN 

42  ANOM  (J)  =CMPLX  (GG  (J)  ,0.) 

CALL  FFTTWO  (ANOM,  NN) 

DO  43  J= 1 , NN 

43  ANOM(J)  =  ANOM  (J)  *B  (J) 

CALL  FFTTWO  (ANOM, NN) 

CALL  REORDR  (A  NOM,  NN) 

DO  44  J=1 , NN 

44  DG  (J) =REAL  (ANOM (J) ) /NN 

IF  (JPRINT  .LT.  3)  WRITE(6,45)  (DG  (J)  ,  J=1 ,  NN) 

45  FORMAT (/////,  30X, 3 1 H  FILTERED  GRAVITATIONAL  ANO M ALY , //, 1 0 ( 


% 
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1  4,/)) 

CC 

cc 

DO  34  KK=1,NKK 
ZUP=Z  Z  (K  K) 

CC 

DO  34  LL=1 , NLL 

IF  (ITOP  .NE.  1)  WRITE  (6,39)  RHO  (LL)  , RHOTWO (LL) 

39  FORMAT (//////, 20X, « BED  THICKNESS  AT  EITHER  OF  PROFILE=»,2F 
DO  40  J= 1 , NN 

40  H  1  ( J)  =0.0 
ISKIP=0 
NTERM  S=3  5 

C  CALL  INVERSION  SBRODTINE  TO  DO  INVERSION 

CALL  DENSTY  (ITOP , RODIG 1 , RHO (LL) , RHOTWO (LL)  ,GG,H1 ,H, 

1SUM, ANOM, HTON, B,SH) 

IF  (JPRINT  .LT.  3)  WRITE  (6r65)  ( RODIG  1  (J)  , J= 1 , N N) 

NTERM  S=  3  5 
ND=NN/2 

C  COMPOTE  GRAVITY  FROM  INVERTED  DENSITY  PROFILE. 

CALL  DENBOU (DG ,RODIG1 ,  Hr SOM , HTON r NTERMS) 

IF  (JPRINT.LT. 3)  WRIT  E  (  6 , 5  1 )  (GG  ( J)  ,  DG  ( J)  ,  J=  1  ,  N  N) 

IF  (JPRINT. GE. 3)  WRITE(6,512)  (GG  (J)  , J=1 , ND) 

512  FORMAT (/////,  30X,  *  ORIGINAL  GRAVITY' ,//, 10 (10F1 2. 

1  *,/)) 

IF  (JPRINT. GE. 3)  WRITE(6,511)  (DG  (J)  ,  J=1 ,  ND) 

511  FORMAT (/////, 30X, '  GRAVITY  FROM  INVERTED  TOPOG ' r // , 1 0 ( 1 0 FI 
1  <+,/)) 

51  FORMAT (/////, 40X, 48H  OBSERVED  AND  CALCULATED  GRAVITATIONAL 

1 1 E  S  , // ,  1 0  (2F9. 3 , 6X  ,  2F9.3,6X,  2F9.3,6X,  2F9.3,6X,  2F9. 

34  CONTINOE 

STOP 
END 

SOBROOTINE  DENSTY (ITOP , RODIG1 , RHOI , RH02,GG, HI , 

1H,SUM, ANOM, HTON, FIL,SH) 

C  ##  *  #  #  #  Jje  *  :)c  #  *  *  #  #  ##  *  *  *  * 

C  ROOTINE  COMPOTES  THE  DENSITY  FROM  RESIDOAL  GRAVITY  PROFILE. 

C  **  ITOP  .EQ.  1  IF  H (X)  AT  EITHER  ENDS  OF  PROFILE  IS  GIVEN ,  ELSE 
C  **  RODIG  1  VECTOR  FOR  DIGITISED  H(X) 

C  **  RHOI  H ( X)  AT  W  END  OF  PROFILE. 

C  **  RH02  H ( X)  AT  E  END  OF  PROFILE. 

C  **  GG  RESIDOAL  GRAVITY. 

C  **  HI  VECTOR  CONTAINS  INITIAL  GUESSED  VALUES  OF  DENSITY  . 

C  **  H  VECTOR  CONTAINS  INTERMEDIATE  AND  FINAL  COMPUTED  VALUES  OF 
C  **  SUM  IS  THE  SUMMATION  OF  THE  FFT  OF  GRAVITY  AND  DENSITY  TERMS 
C  **  HTON  IS  A  DUMMY  COMPLEX  VARIABLE  FOR  THE  DENSITY  D(X),  WHICH 
C  FOR  FFT  OF  DENSITY. 

C  **  ANOM  IS  A  DUMMY  VARIABLE  FOR  CIMPLEXING  THE  RESIDUAL  GRAVITY 
C  USED  FOR  FFT  OF  GRAVITY  DATA. 

C  **  FIL  IS  THE  FILTER  COEFF  SUPPLY  BY  LOPASS  SUBROUTINE. 

C  **  SH  IS  THE  CUT  OFF  FREQUENCY. 

Q  ###****:<£#***  ******  ************  ******  ************************ 

COM MON/HILL S/NN, DX , ZU P , GRH O, G RO 


■ 
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COMMON/INV/NT  ERMS  ,  ITER, ISKIP, ACC ,IZUP,IPRINT,FRAC,NLAY 
C  *  ******************************************  ******************** 

DIMENSION  RODIG1  (1)  ,  Hi  (1)  ,  H  (1)  rGG  (1)  ,  FIL  (1) 

C  **************************************  *************  ********* 
COMPLEX  SOM(1)rHTON  (1)  , A  NOM  (  1 ) 

COMPLEX  HC 

DATA  PI/3.  1415926535/, G/6. 67/ 

CC 

TWO  PI  =  2*PI 
NND2=NN/2 
N2PLUS=NND2+1 
DK=TWOPI/  (DX*NN) 

CC 

CC  CALCULATE  CONTEIBUTIONS  BELOW  THE  CUTOFF  FREQUENCY 
NFINAL  =  NN  *SH  +2 

IF (NFINAL  .GT.  N2PLUS)  NFI NAL=N2 PLUS 

DO  10  J= 1 , NN 

SUM (J) =CMPLX  (0. ,  0.) 

10  ANOM(J)  =CMPLX  (GG(J)  ,0.) 

C 

CALL  FFTTWO(ANOM,NN) 

C 

SOL  D= 1 . E0  5 

C  IF  H (X)  ALONG  PROFILE  IS  NOT  KNOWN  ASSUME  LINEAR  VARIATION  OF 
IF  (ITOP  .EQ.  1)  GOTO  621 
ND2=NN/2 

ROIN  =  ( RHO  1-RH02 )  /  (ND2-  1) 

DO  303  1=1, ND2 

RODIG1  il) =  RHOI -ROIN* (1-1) 

ROD IG 1  (I  +  ND2) =RH01-ROIN*  (ND2-I) 

303  CONTINUE 

Q  ********************************** 

IF  (IPRINT.LE.  2)  WRITE(6,373)  (RODIG  1  (I )  ,1=  1  ,  NN ) 

Q  ************************************  *******  ************ 

621  DO  50  J J= 1 , IT ER 

BIG=- 1 . E09 
SMALL= 1 . E09 
DO  12  J=  1  , NN 

IF(H1  (J)  .GT.  BIG)  BIG=H1  (J) 

IF (HI  (J)  . LT. SMALL)  SMALL=H1(J) 

12  CONTINUE 
WRITE  (6, 13)  JJ, BIG, SMALL 

13  FORMAT (//, 5X, 6H  ITER= , 14, 5 X, 1 3H  MAX  DE NSITY= , F 1 3 . 4 , 5X 

1  , 13  H  MIN  DENSITY  =  , F 1 0 . 4) 

CC  INITIALIZE  SUM 

Q  **************************** 

31  DEXP=EXP  (DK*ZUP) 

TERM=1./(TWOPI*G) 

SUM  (1)  =  ANOM  (1)  *TERM 
DO  15  K=2, NFINAL 
TERM=TERM*DEXP 


. 
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IF  (K.  NE.N2PLUS)  SUM(K)  =  ANOM(K)*TERM 
MINUS=NN-K+2 

15  SUM  (MINUS) =  ANOM  (MINUS)  *TERM 

c  ***************************** 

c 

IF (IPRINT. EQ. 0)  WRITE  (6,55)  (SUM ( J)  , J= 1 , NN) 

55  FORMAT (////, 19H  FOURIER  AMPLITUDES,//,  1 0 ( 1 0 E 1 2 . 3 , /) ) 

C  DIGITISE  DENSITY 

3 73  FORMAT (///,10X, *  DIGITISED  TOPOG=», 

1  //,  1 0 ( 10F 12.4,/) ) 

C  IF  THE  INITIAL  DENSITY  IS  ZERO  THEN  DO  NOT  CALCULATE  THE  PER 
IF  (TSKIP.EQ. 0)  GO  TO  20 
FAC= 1 . 0 

DO  16  N=2 , NTERMS 
FAC=N*FAC 

CC  COPY  THE  NTH  POWER  OF  THE  TOPOGRAPHY  INTO  COMPLEX  HTON 
DO  724  J= 1 , NN 

724  HTON (J) =CMPLX  (HI  (J) *RODIG1 (J) **N, 0. 0) 

127  CALL  FFTTWO (HTON , NN) 

C 

C 

CC  FIND  THE  LARGEST  VALUE  OF  THE  PERTURBATION  TERM 
BIG=0.0 

DO  18  K=2, NFINAL 

CAYN= ( (K-1) *DK) **  (N- 1 ) /FAC 

IF (K. EQ. N2PLUS)  GO  TO  40 

HC=HTON (K) *CAYN 

B=C ABS  (HC) 

IF (B. GT. BIG)  BIG=B 
SUM  (K)  =SUM  (K-HC 

40  MINUS=NN-K+2 

HC= HTON (MINUS) *CAYN 
B=CABS (HC) 

IF (B. GT. BIG)  BIG=B 
18  SUM  (MINUS)  =SUM(MINUS)  -  HC 

IF ( N . EQ . 2)  FIRST=BIG 
IF  (IPRINT.  LE.  2)  WRITE(6,41)  N,  BIG 

41  FORMAT (30X,3H  N  =  ,I4, 1  OX  ,51 H  LARGEST  ABSOLUTE  VALUE  OF  TERM 
1  WARD  PROBLEM=,  E15.4) 

IF (BIG/FIRS T  .LT.  5.E-03)  GO  TO  20 

16  CONTINUE 
WRITE  (6,56) 

56  FORMAT  (20X,34H  THE  FORWARD  SUM  HAS  NOT  CONVERGED) 

C 

CC  APPLY  THE  LOWPASS  FILTER 
20  IF (IPRINT. EQ. 0)  WRITE(6,55)  (SUM  ( J)  , J= 1 , NN) 

Q  J****************  ******  ******** 

C  INVERSE  TRANSFORM  TO  GET  TOPOG,  DIVIDE  BY  CORRESPONDING 

C  TOPOG,  TRANSFORM  AND  FILTER  AND  COMPARE  EACH  ADJACENT 

C  ITERATION. 

CALL  FFTTWO  (SUM, NN) 

CALL  REORDR (SUM, NN) 

DO  272  J= 1 , NN 
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272  H  (J)  =  RE  AL  (SOM  (J)  )  /  (NN*RODTG1  (J)  ) 

DO  299  J= 1 r  NN 

299  SUM(J)  =CMPLX  (H  (J)  ,0.0) 

CALL  FFTTWO  (SOM, NN) 

DO  37  K=1,NN 

3  7  SUM  (K)  =SU M  (K)  *FIL(K) 

CALL  FFTTWO (SUM,NN) 

CALL  REORDR (SUM, NN) 

S1=0.0 

DO  263  1=1, NN 
H  (I)  =REAL  (SUM  (I)  )  /NN 
263  S1=S1+ (H  (I-H1 (I)  )  **2 

C  *************************** ****** 

727  Si =SQRT (SI) /NN 

W RITE  (6,23)  JJ,S1 

23  FORMAT (//, 40  X , 1 1 H  ITERATION  =  , 15 , 1 0X , 7H  ERROR=, El  2 . 4) 

IF  (IPRINT.  LE.  2)  WRITE(6,24)  ZUP,  (H  (  J)  ,  J=  1  ,  NN) 

24  FORMAT (///, 10  X, 4 1 H  THE  DENSITY  MEASURED  RELATIVE  TO  ZUP=, 

1  F  8 . 4 , // ,  10 (10F12. 4,/) ) 

C 

CC  CHECK  TO  SEE  IF  ITERATION  IS  STARTING  TO  DIVERGE 
XX=  (SI-SOLD)  /SOLD 
C 

IF  (XX  .LT.  .01)  GO  TO  42 
WRITE (6 , 43) 

43  FORMAT (///,20X,  57H  ****  THE  ITERATION  PROCEDURE  IS  STARTIN 
1 VERGE  ****,///) 

USE  THE  DENSITY  FROM  THE  LAST  ITERATION 
DO  44  J= 1 , NN 

44  H  (J)  =H  1  (J) 

GO  TO  381 

C  ************************************************* 

42  IF  (SI .LE. ACC)  GO  TO  30 
SOL  D=S 1 

C  ** **  ***  ****  **  * * ***  *****  * ********* 

742  DO  25  J= 1 , NN 

25  HI  (J)  =H  (J) 

50  ISKIP= 1 

C 

381  WRITE (6,26)  S1,ACC,ITER 

26  FORMAT (///,44H  *********  CONVERGENCE  HAS  NOT  BEEN  OBTAINED 

1  1  OX , 7H  ERROR=,E12.3, 10X,18H  DESIRED  ACCURACY  =  , E12. 3, 1  OX, 

1  22H  NUMBER  OF  ITERATIONS=, 15) 

30  IF  (IPRINT.  GT.  2)  WRITE(6,24)  ZUP ,  (H  ( J)  ,  J=  1 ,  NND2) 

RETURN 

END 

SUBROUTINE  DENBOU ( DG, RODIG1 , H, SUM, HTON , NTERMS) 

CC$$$$ $  CALLS  FFTWO 

C  PROGRAM  COMPUTES  GRAVITY  FROM  INVERTED  DENSITY  PROFILE. 

C  DG  CONTAINS  COMPUTED  FINAL  GRAVITY  FROM  INVERTD  DENSITY  PRO 

C  RODIG1  CONTAINS  DIGITISED  H (X)  ALONG  THE  PROFILE. 

C  H  CONTAINS  DIGTISED  DENSITY  VALUES  ALONG  PROFILE. 


' 
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C  SUM  AND  HTON  ARE  THE  SAME  AS  THOSE  OF  DENSTY  SUBROUTINE. 

CC  COMPUTES  2-D  BOUGUER  ANOMALIES  WITH  F FT  POWER  METHOD 
CC  • NTERMS'  IS  THE  NUMBER  OF  TERMS  IN  POWER  SERIES  TO  BE  TAKEN. 
CC  PARAMENTS  COME  THROUGH  COMMON  /HILLS/  -  1 NN '  NUMBER  OF  TERMS 
CC  'DX*  SPACING  IN  X  DOMAIN, «ZUP'  VERTICAL  DISTANCE  ABOVE  Y  =  0  OF 
CC  '  G  R HO '  BIG  G  TIMES  RHO,«H'  ARRAY  WITH  »NN'  ELEMENTS  CONTAININ 
COMMON/HILLS/NN, DX , ZU P, GRH 0, GRO 
DIMENSION  RODIGl  (1)  ,H  (1)  ,DG(1) 

COMPLEX  SUM  (1  ),  HTON  (1)  ,  HC 
COMMON/BOUG/ACC, IPRINT 
DATA  TWOPI/6. 283185307/ 

DATA  G/6.67/ 

WRITE  (6,20)  ACC 

20  FORMAT (///,30X,29H  ACCURACY  LEVEL  IN  BOUGER  IS  =  ,  E  1  2 . 3  , //) 
CC  SETS  VARIOUS  CONSTANTS  AND  ZEROS  1  SUM  *  READY  TO  RECEIVE  SERIE 
NNBY2=NN/2 
N2PLTJS=NNBY2+  1 
DK=TWOPI/(DX*NN) 

DO  110  1=1, NN 

110  SUM  (I)  =(0.0, 0.0) 

FAC=1 .0 
BIGM 1 =10. 

C 

CC  'N'  COUNTS  NUMBER  OF  TERMS  IN  SERIES. 

DO  2500  N=1, NTERMS 
FAC=N*FAC 

CC  COPIES  N-TH  POWER  OF  TOPOGRAPHY  INTO  COMPLEX  • HTON  * 

BIG=0 .0 
DO  233  J= 1 , NN 

233  HTON (J) =C MPLX  (H (J) * RODIGl (J) **N,  0. 0) 

C 

234  CALL  FFTTWO (HTON, NN) 

C 

CC  FILTERS  FOURIER  COEFFS  AND  ADDS  THEM  INTO  'SUM' 

CC  FEWER  THAN  NTERMS  IN  THE  SUM  MAY  BE  TAKEN  IF  AN  ABSOLUTE  CON 
CC  CRITERION  IS  MET. 

CC  ALSO  PERFORM  THE  UPWARD  CONTINUATION  AT  THIS  STAGE 
TERM=TWOPI*G/NN 
DEXP=EXP  (-DK*ZUP) 

IF(N.EQ.I)  SUM  (1)  =HTON  (1)  *TERM 

DO  2200  K=2, N2PLUS 

TERM=TERM*DEXP 

C A Y N=  (  (K-1  )  *DK)  **  (N-1J/FAC 

CAYN=CAYN*TERM 

IF (K. EQ. N2PLUS)  GO  TO  40 

HC=HTON (K) *CAYN 

B=C ABS  (HC) 

IF (B. GT. BIG)  BIG=B 
SUM  (K) =SUM  (K)  +HC 

40  MI NUS  =  NN-K+  2 

HC  =  HTON  (MINUS) *CAYN 
B=CABS (HC) 

IF (B. GT. BIG)  BIG=B 
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2200  SUM (MINUS) =SUM(MINUS)  +HC 
IF(N.EQ.I)  FT  RST=  BIG 
IF  (IPRTNT.LE. 5)  WRITE(6,42)  NrBIG 

42  FORMAT (30X,3H  N  =  , 14, 1  OX  ,  51 H  LARGEST  ABSOLUTE  VALUE  OF  TERM 
1WARD  PROBLEM=,  El  5.  4) 

BB= AM  AX  1  (BIG, BIGM1) 

IF (  BB/FIRST  .LT.  ACC  )  GO  TO  43 
BIGM 1 =BIG 
2500  CONTINUE 

WRITE  (6  , 44)  FIRST, BB 

44  FORMAT (////, 51H  *****  THE  FORWARD  PROBLEM  HAS  NOT  CONVERGE 
1  , 1  OX, 7H  FIRST=, E12.3,10X,5H  BIG=,E12.3) 

43  CALL  FFTTWO (SUM , NN) 

CALL  REORDR  (SUM , NN) 

DO  70  1=1, NN 

70  DG(I)  =  REAL  (SUM(I)  ) 

RETURN 

END 
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